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This is a Semi-Annual Status Report on research conducted between 22 March 1991 and 
21 September 1991 under NASA Grant NAG 5-814, entitled "The Interpretation of Crustal 
Dynamics Data in Terms of Plate Motions and Regional Deformation near Plate Boundaries." 
This grant has supported the research of one Investigator (S. C. Solomon), one Research Staff 
(R. Reilinger), and two Ph. D. students (A. F. Sheehan and C. J. Wolfe) on behalf of the 
NASA Geodynamics and Crustal Dynamics Programs. 

The focus of the research has been in two broad areas during the most recent 6-month 
period: (1) the nature and dynamics of time-dependent deformation and stress along major 
seismic zones, and (2) the nature of long-wavelength oceanic geoid anomalies in terms of 
lateral variations in upper mantle temperature and composition. The principal findings of our 
research to date are described in the accompanying appendices. The first is an abstract of a 
paper given at the SSA Annual Meeting in March 1991, the next three are abstracts of papers to 
be delivered at the Fall AGU Meeting in December 1991, and the last two are reprints of 
recently published papers. 
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APPENDIX 1 


GPS measurements of deformation associated with the 1987 Superstition Hills 
earthquake, Imperial Valley, California: Evidence for conjugate faulting 


by S.C. Larsen, R. E. Reilinger, H. Neugebauer, and W. Strange 


Presented at the 86th Annual Meeting of the Seismological Society of America, 25-27 March 1991, 

San Francisco 

Abstract published in Seismol. Res. Lett., 62, 34, 1991. 
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Large station displacements observed from Imperial Valley GPS campaigns are attributed to 
the November 24, 1987 Superstition Hills earthquake sequence. Thirty sites from a 42 station 
GPS network established in 1986 have been reoccupied during 1988 and/or 1990. 
Displacements at three sites within 3 kilometers of the surface rupture approach 0.5 m. Eight 
additional stations within 20 km of the seismic zone are displaced at least 10 cm. This is the 
first occurrence of a large earthquake (M s 6.6) within a preexisting GPS network. Best-fitting 
uniform slip models of rectangular dislocations in an elastic half-space indicate 130 cm right- 
lateral displacement along the northwest trending Superstition Hills fault and 30 cm left-lateral 
offset along the conjugate northeast trending Elmore Ranch fault. The geodetic moments are 
9.4 x 10 25 dyne-cm and 2.3 x 10 25 dyne-cm for the Superstition Hills and Elmore Ranch 
faults, respectively. Distributed slip solutions using Singular Value Decomposition suggest 
near uniform displacement along the Elmore Ranch fault and concentrated slip to the northwest 
and southeast along the Superstition Hills fault. A significant component of non-seismic 
secular displacement is observed across the Imperial Valley, which is attributed to interseismic 
plate-boundary elastic deformation. 
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APPENDIX 2 


A comparison of STRC90 and STRC91 campaign results: Strain in the 
Coachella Valley, southeastern California 


by Lewis E. Gilbert, Shawn Larsen, Robert Reilinger, John Beavan, Ken Hudnut, Bill Young, and Bill Strange 


To be presented at the Fall Meeting of the American Geophysical Union, San Francisco, 9-13 December 1991 


Abstract published in AGU 1991 Fall Meeting, Eos Trans. AGU, suppl., 1 17, 1991. 
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The consortium of organizations indicated in the author list and symbolized by our 
unpronounceable acronym has been monitoring deformation in the Imperial and Coachella 
Valleys with increasing detail since the establishment of a first epoch in the southern Imperial 
Valley by NGS in 1986. In 1988 the first epoch measurements were repeated and station 
coverage was densified. In addition to reoccupation of the 1988 network, 1989 saw the 
network extend into the Coachella Valley. During 1990, in cooperation with colleagues 
working south of the Mexican border, 103 stations spanning the region from just south of the 
intersection of the San Andreas (SAF) and San Jacinto faults to the head of the Gulf of 
California were occupied. In addition to the primary marks, thirty-one additional marks were 
occupied for short sessions to provide increased station density along the Banning-Mission 
Creek section of the SAF. In 1991, our emphasis was on the Coachella section of the 
network, with few observations in the southern part of the network. 

Comparison of the 1986 and 1988 results shows considerable deformation associated with 
the Superstition Hills earthquake sequence. Displacements at 3 sites within 3 km of the surface 
rupture approach 0.5 m. Eight additional stations within 20 km of the seismic zone are 
displaced at least 10 cm. Geodetic moments for the events are 9.4 x 10 25 dyne-cm and 2.3 x 
10 25 dyne-cm for the Superstition Hills and Elmore Ranch faults, respectively. These values 
are comparable to the calculated seismic moments. 

The 1990 data have been processed and will be compared with measurements in the south 
and with the 1991 data from the Coachella Valley. 
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APPENDIX 3 


Lateral variations in upper mantle temperature and composition inferred from 
shear-wave travel time and attenuation, geoid, and bathymetry 


by Anne F. Sheehan and Sean C. Solomon 

To be presented at the Fall Meeting of the American Geophysical Union, San Francisco, 9-13 December 1991 
Abstract published in AGU 1991 Fall Meeting , Eos Trans . AGU> suppL, 514, 1991, 
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By means of a large data base of digital seismograms and waveform cross-correlation and 
spectral ratio techniques, we have measured SS-S differential travel time residuals and 
differential attenuation in order to determine lateral variations in upper mantle structure beneath 
the Mid-Atlantic Ridge and East Pacific Rise. After removing the signature of lithosphere age, 
we find evidence for variations in travel time residuals and attenuation along both ridge systems 
at wavelengths of 1000 to 7000 km. These travel time anomalies correlate qualitatively with 
along-axis variations in bathymetry, geoid height, and Love wave phase velocity. We 
formulate a joint inversion of SS-S travel time residual, geoid height, and bathymetry under 
that assumption that all arise from variations in upper mantle temperature or bulk composition 
(parameterized in terms of Mg#). Analysis of differential attenuation places further constraints 
on the lateral heterogeneity of the temperature field. Inversion for temperature perturbations 
alone provides good fits to travel time and geoid variations in the north Atlantic region. 
Compositional variations alone are unable to match the travel time and geoid or bathymetry data 
simultaneously. Temperature variations are ± 50 K and compositional variations are ± 0.5-3 % 
Mg# for models with the temperature variations uniformly distributed over the uppermost 300 
km and the compositional variations either distributed uniformly over the same interval or 
concentrated at shallower depths. The magnitudes of these variations along the Mid-Atlantic 
Ridge are consistent with the chemistry and geothermometry of dredged peridotites and with 
variations in the depth extent of mantle melting implied by the chemistry of dredged basalts. 
Differential travel times of SS-S pairs in the vicinity of the East Pacific Rise show less age 
dependence than residuals in the north Atlantic and are more consistent with the presence of 
azimuthal anisotropy. These contrasts may arise from differences in the dynamics of slow and 
fast spreading ridges. 
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APPENDIX 4 


Oceanic transform earthquakes with unusual mechanisms or locations: Relation to 
fault geometry and state of stress in the lithosphere 


by Cecily J. Wolfe, Eric A. Bergman, and Sean C. Solomon 


To be presented at the Fall Meeting of the American Geophysical Union, San Francisco, 9-13 December 1991 


Abstract published in AGU 1991 Fall Meeting, Eos Trans. AGU, suppl., 518, 1991. 
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On oceanic transforms, most earthquakes are expected to occur on the principal transform 
displacement zone (PTDZ) and to have strike-slip mechanisms consistent with offset of a 
spreading ridge. We conducted a search along the global ridge system for transforms departing 
from this pattern using the Harvard centroid moment tensor catalogue and earthquake locations 
from the International Seismological Centre. Events with unusual mechanisms occur on 
several oceanic transforms. We have analyzed source mechanisms and centroid depths of 
earthquakes on the St. Paul’s, Marathon, Owen, Heezen, Tharp, and Menard transforms from 
an inversion of long-period body waveforms. Relative locations of earthquakes along these 
transforms have been determined with a multiple-event relocation technique. 

Much of the anomalous earthquake activity can be associated with the presence of complex 
fault geometry or large structural features that influence the fault. Reverse-faulting earthquakes 
occur at the St. Paul’s transform near St. Paul’s Rocks and at a compressional bend at the 
Owen transform in the area of Error Seamount. An extensional event on the Heezen transform 
is located at the edge of an extensional offset of the transform, and extensional activity along 
the Tharp and Menard transforms may also be correlated with extensional offsets. In general, 
prominent gaps in seismicity can exist near both extensional and compressional bends in a 
transform. 

Some earthquakes with unusual mechanisms occur outside of the transform fault zone and 
do not appear related to fault zone geometry. Earthquakes with thrust mechanisms are located 
beneath transverse ridges near the ridge-transffom intersection of the St. Paul’s transform and 
the Marathon transform. A normal-faulting earthquake occurred outside the Heezen transform 
fault zone, in an elongate transform-parallel basin. 

Possible influences on the occurrence of such earthquakes include lithospheric thermal 
stresses, tractions at the base of the plate from asthenospheric flow, changes in plate motion, 
wrench tectonics, inherited structure, and rotation of principal stresses at a weak fault. 
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Age constraints for the present fault configuration in the Imperial Valley, 
California: Evidence for northwestward propagation of the Gulf of California rift 

system 


by Shawn Larsen and Robert Reilinger 
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Age Constraints for the Present Fault Configuration in the Imperial Valley, 
California: Evidence for Northwestward Propagation of the 
Gulf of California Rift System 

Shawn Larsen 

Seismological Laboratory, California Institute of Technology , Pasadena 

Robert Reilinger 

Earth Resources Laboratory, Massachusetts Institute of Technology , Cambridge 


Releveling and other geophysical data for the Imperial Valley of southern California suggest 
the northern section of the Imperial- Brawley fault system, which includes the Mesquite Basm and 
Brawley Seismic Zone, is much younger than the age of the valley itself. A minimum age of 3000 
years is calculated for the northern segment of the Imperial fault from correlations between surface 
topography and geodetically observed seismic/interseismic vertical movements. Calculation of a 
maximum age of 100,000 years is based upon displacements in the crystalline basement along 
the Imperial fault, inferred from seismic refraction surveys. This young age supports recent 
interpretations of heat flow measurements and the evolution of geothermal systems, which also 
suggest that the current patterns of seismicity and faulting in the Imperial Valley are not long 
lived. The current fault geometry and basement morphology suggest a northwestward growth of 
the Imperial fault and a northwestward migration of the Brawley Seismic Zone. If this localized 
process is representative of more regional tectonic processes along the extent of the Salton Trough, 
we suggest that this migration is a manifestation of the propagation of the Gulf of California rift 
system into the North American continent. 


Introduction 

The Salton Trough is a complex transition zone between 
crustal spreading in the Gulf of California and right-lateral 
transform motion along the San Andreas fault system (Fig- 
ure 1). The Imperial Valley is that section of the Salton 
Trough north of the U.S. -Mexico border and south of the 
Salton Sea (Figure 2). The Trough is characterized by pre- 
dominately right-stepping, right-lateral en echelon faults, 
presumably linked by zones of crustal extension [Lomnitz et 
al, 1970; Elders et of., 1972]. It forms a 150 by 300 km struc- 
tural depression which is filled by up to 15 km of late Ceno- 
zoic sediments. The seismic velocity of the lower 5-10 km 
(V p = 5.7 km/s) suggests these sediments are greenschist- 
facies, metasedimentary rocks [Fuis et al., 1984]. The age of 
the Imperial Valley-Salton Trough region is suggested to be 
between 4 and 12 million years [Larson et al., 1968; Moore 
and Buffington, 1968; Ingle, 1974]. 

The Imperial Valley and its major fault systems trend 
northwesterly, nearly parallel to the relative motion between 
the North American and Pacific plates. Dextral faulting pre- 
dominates, although northeast-trending sinistral structures, 
as well as dip- slip motion along north-south surface breaks, 
play a significant role in the regional tectonics [Johnson and 
Hutton , 1982; Nicholson et al, 1986; Reilinger and Larsen , 
1986]. 

The Mesquite Basin is a subaerial topographic low 
bounded on the west by the northern Imperial fault and on 
the east by the Brawley fault (Figure 2). Maximum basin 
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relief is about 10 m relative to its periphery. Evidence that 
the Mesquite Basin is actively subsiding includes geodetic 
measurements of surface deformation and measurements of 
vertical slip along the Imperial and Brawley faults. We pro- 
vide evidence that the Mesquite Basin is extremely young 
compared to the age of the Imperial Valley, suggesting this 
section of the Imperial-Brawley fault system is at an early 
stage of tectonic development. We extend this hypothesis 
and suggest ongoing northwestward propagation of the Gulf 
of California rift system. 

Imperial Valley Seismicity and Faulting 

The Imperial Valley is one of the most seismically active 
regions of California (Figure 3). A significant fraction of 
this seismicity occurs within the Brawley Seismic Zone, a 
region of high activity between the northern Imperial and 
southern San Andreas faults [ Johnson , 1979; Johnson and 
Hill , 1982]. The Imperial fault ruptured historically in 1940 
(Ml 6.4, Ms 7.1) and in 1979 (Ml 6.6, M s 6.9); episodes of 
creep have been recognized along the fault since 1966 [Allen 
et al, 1972]. The seismic moment of the 1979 earthquake 
is well determined (6 x 10 25 dyn cm) [e.g., Kanamori and 
Regan , 1982], while that for the 1940 event ranges between 
10 and 80 x 10 25 dyn cm [Trifunac and Brune , 1970; Hanks 
et al, 1975] although the 48 X 10 25 dyn cm moment esti- 
mated from long-period surface waves is preferred [Doser 
and Kanamori, 1987]. Other major earthquakes in the 
Imperial Valley include the recent 1987 Superstition Hills 
earthquake sequence: a Ms 6.2 event produced by slip along 
a northeast-trending seismic lineament, followed 12 hours by 
a Ms 6.6 earthquake produced by slip along the Superstition 
Hills fault [Magistrate et al., 1989; Williams and Magistrate, 
1989]. 

The 1979 surfacial rupture of the Imperial fault extended 
from a point 5 km north of the border northwestward 
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transform motion along the San Andreas fault. The Imperial Val- 
ley is that portion of the Salton Trough north of the U.S. -Mexico 
border and south of the Salton Sea. Abbreviations are S.F., San 
Francisco; L.A., Los Angeles. Map modified from Lachenbruch 
et al. [1985]. 


33.1 km to a point south of Brawley (Figure 4). The pre- 
dominate strike of the Imperial fault is N37° W. Along 
the northwestern most 5 km, however, the fault bends and 
trends north; we refer to this segment as the north exten- 
sion. Trending parallel and lying 6 km east of the north 
extension, the Brawley fault ruptured in 1979 along a 13 km 
surface break. This rupture pattern generally featured left- 
stepping en echelon cracks that extended a few millimeters 
to a few centimeters [Sharp et al., 1982]. A third, relatively 
minor 1 km north-trending break named the Rico fault was 
mapped 6-7 km east of the Brawley fault (Figure 4). The 
surface breakage along this structure resembled that of the 
Brawley fault zone. The geometrical similarity in strike and 
separation shown by the north extension, Brawley, and Rico 
faults, suggest a similar tectonic origin. 

The epicenter of the 1940 earthquake was north of the 
U.S. border, but right-lateral surfacial offsets were larger 
in Mexico (Figure 3). A maximum surface offset of 6 me- 
ters was recorded near the border, with displacement ta- 
pering off rapidly to the north [Trifunac and Brune, 1970; 
Sharp , 1982]. Geodetic measurements indicate 4.5 and 3.0 m 
of right-lateral slip (coseismic plus postseismic) along the 
southern and northern halves of the Imperial fault, respec- 
tively (i.e., north and south of the epicenter), with 2.0 m 
postseismic slip along a northwest extension of the Brawley 
fault [ Reilinger , 1984]. The 1979 epicenter was south of the 
border, although surfacial displacement was observed only in 
the United States. Maximum coseismic surfacial offset was 
55-60 cm, with considerable afterslip (~ 30 cm) during the 
following 6 months [Sharp et al., 1982]. Strong ground mo- 
tion and geodetic modeling [ Archuleta , 1984; Hartzell and 
Heaton, 1983; Reilinger and Larsen , 1986] suggest an av- 
erage slip of about 1 m along the fault plane, with small 
patches of greater displacement (inferred to be asperities). 

The mechanism of strain transfer between the Imperial 
and San Andreas faults within the Brawley Seismic Zone has 
been the focus of considerable investigation [e.g., Johnson , 
1979]. A conjugate relationship of right-lateral, north west- 



Fig. 2. The Imperial Valley and important faults. The Mesquite Basin is a subaerial topographic depression of 
about 10 m between the Imperial and Brawley faults. 
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SEISMICITY 1932-1989 
CALTECH CATALOG 



RELOCATED 1979 
AFTERSHOCKS 



Fig. 3. Seismicity in the Imperial Valley between 1932 and 1989 (Caltech/USGS Catalog) Major events mclude 
the 1940 Imperial Valley (M s 7.1), 1968 Borrego Mountain (M L 6.5), 1979 Imperial Valley (M s 6.6), and the 
1987 Superstition Hills (Ms 6.6, Ms 6.2) earthquakes. The Brawley Seismic Zone is the active region between 
the northern reach of the Imperial fault and the southern extent of the San Andreas. The Sand Hills Seismic 
Lineament is shown by the shaded strip outlining earthquakes trending southeast from the southern end of the 
San Andreas fault. Shown in the inset are aftershocks of the 1979 earthquake which have been relocated following 
the methods outlined by Doser and Kanamor, [1986], The dashed lines represent orthogonal faults used to satisfy 
the observed vertical deformation from the 1979 event [Reilinger and Larsen, 1986]. Focal mechanisms (lower 
hemisphere, equal area projections [Reasenberg and Op P enhtimer, 1985]) for events defining a northwest trend 
indicate right-lateral strike slip motion. 



Fig. 4. Map of the Imperial Valley and important tectonic fea- 
tures. Abbreviations are RF, Rico fault; BF, Brawley fault; NE, 
North Extension. The shaded pattern along each fault indicates 
the surface rupture from the 1979 earthquake. The Brawley Seis- 
mic Zone (hatched) is the region of high seismicity extending 
northwest from the northern reach of the Imperial fault. The 
Salton Sea and Brawley geothermal fields are indicated by the 
shaded patterns. Refraction surveys [Fuis et al., 1984] cross the 
Imperial fault at RL-1, RL-2, RL-3. The leveling route is shown 
by the series of dots from the central Imperial Valley to the east- 
ern border of the Salton Sea (each dot representing a benchmark). 


trending faults perpendicular to left-lateral, northeast 
trending structures may play a significant role in the re- 
gional tectonics [Nicholson et al., 1986]. Although the Impe- 
rial and San Andreas faults strike predominately northwest 
(right-lateral), a left-lateral structure extending northeast 
from the northern terminus of the Imperial fault is indi- 
cated from focal mechanisms and the aftershock pattern of 
the 1979 earthquake [Johnson and Hutton , 1982]. A conju- 
gate fault mechanism is supported by Reilinger and Larsen 
[1986], who suggested several tectonic models of the Brawley 
Seismic Zone satisfying geodetically determined measure- 
ments of vertical surface displacement. The preferred model 
consists of a northeast-trending left-lateral fault conjugate 
to a right-lateral northwest-trending fault dipping 70° to the 
southwest (Figure 3, dashed lines). Neither fault broke the 
surface, but roughly 1 m of slip at depth is required to fit the 
geodetic measurements. A similar conjugate fault relation- 
ship was observed for the 1987 Superstition Hills earthquake 
sequence [e.g., Magistrate et al, , 1989]. 

Aftershocks from the 1979 earthquake have been relo- 
cated following the methods of Doser and Kanamori [1986] 
and Klein [1985] (Figure 3). The northeast-trending seismic 
lineament first identified by Johnson and Hutton [1982] is 
clearly defined. To the north, a tightly constrained group 
of events following a northwesterly direction is indicated. 
Epicentral depths for this cluster range from 5 to 11 km, 
possibly putting them on the 70° west-dipping structure 
suggested by Reilinger and Larsen [1986]. We have com- 
puted focal mechanisms for these events and find them con- 
sistent with a northwest-trending right-lateral fault (Fig- 
ure 3). Thus, both seismic and geodetic data suggest the 
tectonic framework of the Brawley Seismic Zone is marked 
by en echelon northwest-trending right-lateral faults linked 
by conjugate left-lateral structures. 

Extending southeast from the southern tip of the San An- 
dreas fault is a linear alignment of earthquakes [e.g., Johnson 
and Hutton, 1982], here referred to as the Sand Hills Seis- 
micity Lineament (Figure 3). This feature may signify the 
southeasterly extension of the San Andreas fault. Although 
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there is no surfacial geological evidence to support this hy- 
pothesis [Sharp, 1982], creep is suggested along northwest 
segments of the Sand Hills lineament [Jennings, 1975]. 

The earthquake recurrence interval along the Imperial 
fault is not well constrained. Sykes and Nishenko [1984] use 
the 39 year interval separating the 1940 and 1979 shocks as 
well as a 1915 earthquake sequence located near El Centro 
[Beal, 1915] to estimate a 32 year recurrence rate. Anderson 
and Bodin [1987] suggest the fault north of the border will 
next rupture between 2010 and 2050 (50 year recurrence), 
and the next break along the southern segment to occur be- 
tween 2170 and 2290 (300 year recurrence). Measurements 
of surface offset, as well as seismic and geode tically deter- 
mined estimates of fault slip at depth, indicate the 1940 fault 
rupture was several times larger than in 1979, in agreement 
with the larger moment for the 1940 event. North of the 
border, however, the magnitude of horizontal surface dis- 
placement was relatively constant for the two earthquakes. 
One explanation for this is that the fault north of the border 
may rupture more frequently but with smaller events. Alter- 
natively, the large postseismic slip following the 1940 earth- 
quake indicated by geodetic data, suggests that a significant 
fraction of strain buildup may be relieved aseismically. 

If the entire 49 mm/yr movement between the Pacific and 
North American plates predicted by new global plate mod- 
els (NUVEL-1) [DeMets et ai, 1987, 1990] is accommodated 
across the Imperial fault, 1.0 m of seismic or aseismic fault 
slip would require a 20 year interval of strain buildup. More 
likely, however, a significant component of plate motion is 
distributed along the Elsinore and San Jacinto fault systems 
[Sharp, 1981; Pinault and Rockwell, 1984; Snay et at, 1986], 
as well as faults off the coast of southern California [e.g., 
Weldon and Humphreys, 1986], Trilateration and triangula- 
tion measurements from 1941 to 1987 in the central Imperial 
Valley indicate an average displacement across the Imperial 
fault of 35-43 mm/yr [Prescott et al., 1987; Snay and Drew, 
1988]. Preliminary results utilizing the Global Positioning 
System (GPS) suggest a slightly larger rate between 1986 
and 1989, although interpretation of these measurements 
have been complicated by large displacements from the 1987 
Superstition Hills earthquake sequence [Larsen, 1991]. 

Assuming 40 mm/yr of plate motion across the Imperial 
fault, 1.0 m of potential slip will accumulate in 25 years. 
This is equivalent to the earthquake recurrence interval, at 
least for the northern segment of the Imperial fault, if the 
^ 1.0 m surface displacement measured in 1940 and 1979 
is characteristic of fault displacement and if all slip is gen- 
erated seismically. Considering the likelihood of aseismic 
deformation, as well as seismic and geodetic models indicat- 
ing 2-3 m slip asperities along the 1979 rupture plane, it is 
reasonable to expect that the average slip generated along 
the northern Imperial fault during each earthquake cycle 
is somewhat greater than 1.0 m. Assuming 2-3 m of slip 
(based on the seismic plus postseismic offset estimated for 
the 1940 earthquake and the maximum slip observed for the 
1979 earthquake), more reasonable estimates of earthquake 
recurrence would be 50 to 75 years for this segment of the 
Imperial fault. 

Subsidence of the Mesquite Basin 

First-order leveling surveys crossing the northern 
Mesquite Basin were conducted by the National Geodetic 
Survey (NGS) in 1931, 1941, 1974, 1978, and 1980 (Fig- 
ure 4). Profiles of elevation change from 1931 to 1941, 1941 


to 1974, and 1978 to 1980 are shown in Figure 5. The pro- 
cedure used to determine these crustal movement profiles 
is described in Brown and Oliver [1976]. Briefly, an esti- 
mate of relative elevation change between successive bench- 
marks is obtained by subtracting the elevation difference 
between benchmarks measured at some reference time from 
the difference measured at some later time. These movement 
profiles have not been connected to any external reference. 
Therefore, only relative movements along the level lines are 
significant. 


a) Mesquite 

Basin 

I 1 



Fig. 5. (a) Elevation (dashed line) along the releveling route 

between El Centro and the Salton Sea. The adjusted topography 
(solid line) is the elevation with the northward tilt of —0.0011 
radians removed. The 10m depression between 9 and 22 km is the 
surfacial expression of the Mesquite Basin, (b) Elevation changes 
along the leveling route from 1931 to 1941, 1941 to 1974, and 
1978 to 1980. Note the strong correlation between deformation 
and the surface expression of the Mesquite Basin. 

The random error for these measurements is compara- 
tively small, less than 1 cm. In addition, elevation-correlated 
errors (i.e., rod calibration and atmospheric refraction) 
which can obscure or be mistaken for real tectonic deforma- 
tion, will not seriously affect the data because of negligible 
topographic variation along the leveling route (Figure 5a). 

The 1931-1941 and 1941-1974 profiles have been mod- 
eled as coseismic and postseismic deformation from the 
1940 Imperial Valley earthquake [Reilinger, 1984]. Displace- 
ments for the most recent interval (1978 to 1980) have been 
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modeled as surface deformation from the 1979 earthquake 
[Reilinger and Larsen , 1986]. The most striking feature of 
the leveling data is the similar pattern of subsidence across 
the Mesquite Basin observed on all three profiles, suggesting 
this deformation style is characteristic for the region. Co- 
seismic subsidence for the 1940 and 1979 events are on the 
order of 10-15 cm, with an additional 15 cm following the 
1940 earthquake. Total subsidence for the period 1931 to 
1980 is about 40 cm. 

Elevation along the leveling route is shown in Figure 5a 
(dashed line). A relatively constant northward slope of 
—0.0011 radians is observed. This long-wavelength trend 
may mask small scale variations, so we construct a modified 
topographic profile by removing this regional slope (we add 
0.0011 radians to the true profile). The modified profile, or 
adjusted topography, is shown as the solid line in Figure 5a. 
The 10-meter depression between 9 and 22 km marks the 
boundary and surface relief along the northern part of the 
Mesquite Basin. The topographic relief is well correlated 
with the seismically generated subsidence, strongly suggest- 
ing the Mesquite Basin formed by many episodes of seismic 
activity similar to the 1940 and 1979 events. 

Vertical surface slip along the northern section of the 1979 
rupture plane ranged from 0 to 30 cm (including 6 months 
afterslip), while vertical offset along the Brawley fault was 
0 to 24 cm [Sharp et al , 1982]. Measurements of vertical 
slip following the 1940 earthquake were sparse, although the 
sense of displacement was generally the same as in 1979 
[Sharp, 1982]. During an earthquake swarm in 1975, up 
to 20 cm of vertical displacement was observed along the 
Brawley fault and an additional 20 cm possibly occurred 
between 1960 and 1975 [Sharp, 1976]. In each case, slip was 
down to the east along the Imperial fault and down to the 
west along the Brawley fault. Displacement on the Rico 
fault during the 1979 event was down to the west. 

Perhaps the most puzzling and intriguing aspect of de- 
formation in the Mesquite Basin is shown by the offset pat- 
tern recorded in the crystalline basement along the Impe- 
rial fault. Seismic refraction experiments were conducted 
by the U.S. Geological Survey in the Imperial Valley dur- 
ing 1979 [Fuis et al , 1984]. Three refraction lines RL-1, 
RL-2, and RL-3 cross the Imperial fault where shown in 
Figure 4. (These correspond to Fuis et al [1984] lines 6NW- 
1SE-1NW, 1ESE, 1E-2W.) The seismic measurements in- 
dicate a 1000 m basement offset across the Imperial fault 
at RL-1, a 500 m offset at RL-2, whereas no basement off- 
set is observed at RL-3. That is, the offset increases to 
the southeast. Where detectable, the subsurface morphol- 
ogy is down to the east. The basement is defined as rock 
with V P = 5.6 km/sec, which approximately corresponds to 
5 km depth. What makes the basement structure so un- 
usual is its opposite arrangement to the deformation dis- 
played at the surface, where vertical fault offsets measured 
for the 1940 and 1979 earthquakes generally increased to the 
northwest. In fact, where the basement structure is maxi- 
mum (at RL-1), the coseismic vertical surface displacements 
were either small or non-existent. Presumably, this appar- 
ent discrepancy between surface and sub-surface structure 
must illustrate an important tectonic feature. 

Age of Faulting 

The correlation between geodetically measured subsi- 
dence and the topographic expression shown in Figure 5 
strongly suggests this region developed from episodes of seis- 


mic activity similar to the 1940 and 1979 earthquakes. In 
fact, this example clearly illustrates that earthquakes are 
a fundamental building block of tectonic structures. The 
10 m surface depression, together with the subsidence rate 
and basement morphology, places constraints on the age of 
the Mesquite Basin, and correspondingly the northern seg- 
ment of the Imperial fault. 

About 5 m of seismic and postseismic slip along the Im- 
perial fault north of the border is required to form the 
40 cm subsidence between the earliest and most recent level- 
ings across the Mesquite Basin (1931-1980) [Reilinger, 1984; 
Reilinger and Larsen , 1986]. At a slip rate of 40 mm/yr 
across the Imperial fault, 5 m of potential slip will accu- 
mulate in 125 years. The equivalent basin subsidence rate 
is thus about 3 mm/yr. While depending heavily on the 
rate of strain accumulation, this analysis is invariant to the 
earthquake recurrence interval. 

At a tectonic subsidence rate of 3 mm/yr, the 10 m de- 
pression which outlines the Mesquite Basin would form in 
3000 years. This suggests that the tectonic framework un- 
derlying the basin, namely the northern Imperial and Braw- 
ley faults, is extremely young compared to the 4 to 12 million 
year age of the Imperial Valley, However, this estimate does 
not include sediment influx into the Mesquite Basin. While 
the measured seismic subsidence is about an order of mag- 
nitude larger than typical fill rates in arid regions [Ollier, 
1981], the basin is located in one of the largest river deltas 
in the United States; presumably sediment influx is high. In 
fact, the average rate of deposit in the central Imperial Val- 
ley is about 1 mm/yr (5 km over the last 5 million years), 
only slightly smaller than the rate of tectonic subsidence. 
Overlying sediments may mask a deeper basin, so 3000 years 
is an extreme minimum duration for basin development. 

The lack of an observed basement offset at RL-3 places 
further constraint on fault age. The geometry and dextral 
motion of the San Andreas and Imperial faults require ex- 
tension in the Brawley Seismic Zone. Dip-slip motion along 
the northern Imperial and Brawley faults helps to fill this re- 
quirement. Although geodetic, geologic, and strong-motion 
data indicate significant vertical displacements along the 
northern segment of the Imperial fault (north of its inter- 
section with the Brawley fault), apparently insufficient time 
has elapsed to allow the formation of a detectable basement 
offset at its northern extent. The lack of offset suggests this 
region formed relatively recently and is at its earliest stage 
of tectonic development. Fuis et al [1984] suggest that on 
intersecting refraction lines in the Imperial Valley, structural 
boundaries to about 5 km depth agree to within a few tenths 
of a kilometer. Assuming the refraction data can resolve off- 
sets of 300 m (about 1/2 of the offset measured at RL-2), 
at a tectonic subsidence rate of 3 mm/ yr the maximum age 
for the northern Imperial fault is about 100,000 years; again 
very young compared to the 4-12 million year age of the 
Imperial Valley. 

Other evidence support a young age for this segment of 
the Imperial fault. Models of heat transfer mechanisms 
suggest the Sal ton Sea geothermal field (Figure 4) formed 
within the last 3000 to 20,000 years [Kasameyer et al, 1980, 
1984], consistent with the 3000 to 100,000 year age range cal- 
culated for the Mesquite Basin. If representative of central 
Imperial Valley tectonics, this geothermal field likely formed 
contemporaneously with the Brawley Seismic Zone and the 
northern Imperial fault. To achieve a balance between ther- 
mal constraints and the current composition of the crust, 
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Lachenbruch et al. [1985] calculate that heat flow measure- 
ments within the Imperial Valley indicate an average exten- 
sion rate of ~ 10~ 14 s" 1 since the formation of the Salton 
Trough. At this rate, the differential velocity between the 
Pacific and North American plates requires that extension 
and faulting must have been distributed over a relatively 
wide region (~ 150 km) during the last several million years. 
Presumably, tectonic and seismic activity, which is presently 
highly concentrated along the Imperial fault and within the 
Brawley Seismic Zone, is part of an evolutionary process in 
which tectonic activity is shifted from one region of the val- 
ley to another. The northern Imperial and Brawley faults, 
Mesquite Basin, and Brawley Seismic Zone may represent 
the most recent epoch of activity in a rapidly changing fault 
geometry. 

Propagating Rift 9 

As discussed above, the relationship between seismicity, 
dip-slip faulting, and basement offset indicates a young age 
for the northern segment of the Imperial fault. Similarly, 
the large basement offset along the central section of the 
fault (at RL-1) suggests significant vertical slip along this 
segment in the past. We suggest here a scenario for the 
recent history of the Imperial fault and the Brawley Seismic 
Zone, that includes the northwestward propagation of the 
Gulf of California oceanic rift system. 

Although rupture along the Imperial fault is predomi- 
nately strike slip, the large component of normal motion 
along the northern segment of the fault is presumably in re- 
sponse to the en echelon geometry of the San Andreas and 
Imperial faults. These faults may act as transforms associ- 
ated with a spreading center beneath the Brawley Seismic 
Zone [Elders et ai, 1972; Johnson, 1979]. If the northern 
extent of the Imperial fault, as well as the Brawley Seismic 
Zone were previously further south (perhaps southeast of El 
Centro), dip-slip motion would be expected along this seg- 
ment of the fault. Eventually, a detectable offset would de- 
velop in the crystalline basement. As the spreading center 
(Brawley Seismic Zone) migrated northwest to its present 
position, so would the vertical movements during seismic 
events. Although rupture along the fault becomes increas- 
ingly strike slip with age, the vertical offset equals the in- 
tegrated offset through time, and therefore increases with 
age. This model accounts for the apparent disparity between 
long-term vertical offsets on the Imperial fault (increasing 
basement offset to the southeast) and present-day seismic 
fault slip (maximum dip-slip along the northern segment of 
the fault). 

The rupture pattern for the 1979 earthquake supports this 
hypothesis (Figure 4). Clearly the northern Imperial and 
Brawley faults are active components in the stress/strain 
transfer mechanism between the Imperial and San Andreas 
faults. Both structures show significant seismic displace- 
ments at the surface and at depth. Although displacement 
along the Rico fault was 10-20 cm vertical with no horizontal 
offset [Sharp et a/., 1982; Reilinger and Larsen, 1986], the 
1 km rupture length suggests it is only a minor constituent 
in the regional tectonics. The Rico, Brawley, and the north 
extension of the Imperial fault, each follow a north-south 
trend and are uniformly spaced at distances of 6 to 7 km. 
This geometrical similarity suggests a similar tectonic origin. 
In fact, the Rico fault may be an older structure reactivated 
during the 1979 earthquake, although little is known about 
its past history [Sharp et al. t 1982]. 


A schematic illustration of the temporal evolution of this 
region is shown in Figure 6. If the Brawley Seismic Zone 
was further southeast than at present, the Rico fault may 
have acted as the Brawley fault does today. Similarly, the 
Brawley fault would have been the northern splay of the 
Imperial fault, identical to the present north extension. A 
prehistoric basin would have developed between the Rico 
and Brawley faults (forming the observed fault offset), simi- 
lar to the Mesquite Basin. Presumably, as the Imperial fault 
lengthened to the northwest, the Rico- Brawley fault system 
no longer influenced the stress/strain distribution between 
the northern Imperial and southern San Andreas faults. As 
a result, a new fault developed (north extension) and the 




Fig. 6. Schematic diagram of past and present fault configurations 
in the Imperial Valley illustrating the hypothesized northwesterly 
migration of the Brawley Seismic Zone. Li this model the Sand 
Hills Seismicity Lineament is the extension of the San Andreas, 
left dormant after the passage of the Brawley Seismic Zone. 
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Rico fault died out. Continued migration of the Brawley 
Seismic Zone may in the future create a new north-south 
trending structure northwest of the present terminus of the 
Imperial fault. As the Brawley Seismic Zone shifted to the 
northwest, so did the southern terminus of the San Andreas 
fault (Figure 6). The Sand Hills lineament appears to be 
the remnant of an older segment of the San Andreas, and 
except for residual seismic activity, left dormant with the 
northwest passage of the Brawley Seismic Zone. 

An apparent inconsistency with this model is that the slip 
direction on the north extension is down to the east, while 
slip along the Brawley and Rico faults is down to the west. 

A distinction is made here between the northern Imperial 
fault and the north extension. It is the normal movement 
on the Imperial fault which generates the offset observed 
in the crystalline basement; the north-trending faults are 
subsidiary features. Once the Brawley Seismic Zone mi- 
grates past its present position, the north extension could 
maintain its current slip orientation without affecting the 
tectonic model suggested in Figure 6. It is also conceivable 
that the slip orientation on this fault could be reversed in 
the near future (down to the west), or that another fault 
could develop in its place. 

It is possible to make a rough estimate for the migra- 
tion rate of the Imperial fault and southern Brawley Seismic 
Zone. Assuming a dip-slip offset rate of 3 mm/yr (estimated 
above), approximately 330,000 years are required to create 
the 1000 m basement offset measured along the Imperial 
fault at RL-1. The 3000 to 100,000 year age for the fault 
segment 20 km to the northwest (at RL-3) indicates that 
the Brawley Seismic Zone has migrated about 20 km during 
the last 250,000 to 300,000 years. This yields a migration 
rate of about 7 cm/yr. While this rate is only a very crude 
estimate, it is significant to note that it is comparable to 
typical plate motion velocities (i.e., several centimeters per 
year). In fact, the estimated spreading rate in the mouth of 
the Gulf of California averaged over the last 3 million years 
is 4.9 cm/yr [e.g., DeMets et al., 1987]. 

If the Brawley Seismic Zone represents the crustal man- 
ifestation of a subcrust al spreading center, we speculate 
that its northwesterly migration is directly associated with 
the propagation of the Gulf of California rift system into 
the North American continent. This hypothesis assumes 
that the localized transient phenomenon observed along the 
northern Imperial fault is characteristic of more regional tec- 
tonic processes throughout the Sal ton Trough. Clearly, the 
Imperial Valley and Brawley Seismic Zone are undergoing 
rapidly changing tectonics. Understanding the kinematics 
of these changes will help constrain the dynamic processes 
which control the transition between crustal spreading in 
the Gulf of California and transform motion along the San 
Andreas fault. 

Conclusions 

Geodetic, seismic, tectonic, and heat flow data in the 
Imperial Valley suggest that the northern segment of the 
Imperial-Brawley fault system is extremely young compared 
to the 4 to 12 million year age of the Imperial Valley. We find 
a minimum age of 3000 years based upon the relationship 
between topography and earthquake induced geodetic dis- 
placements, and estimate a maximum age of 100,000 years 
based upon observed basement offsets across the Imperial 
fault determined from seismic refraction surveys. A young 
age is consistent with heat flow data, which indicate a dis- 


tributed and ephemeral pattern of faulting in the Imperial 
Valley [Lachenbruch et al., 1985]. 

In addition, we speculate that the apparent incompati- 
bility along the Imperial fault between the observed seismic 
vertical displacements (maximum to the north) and the off- 
set recorded in the crystalline basement (maximum to the 
south) is a direct result of the northwestward propagation 
of the Imperial fault and Brawley Seismic Zone. A series 
of evenly spaced north-trending surface ruptures and the 
Sand Hills seismicity lineament are consistent with this hy- 
pothesis. A 7 cm/yr migration rate is calculated from mea- 
sured surface displacements and from variations in basement 
morphology along the Imperial fault. The migration of the 
Brawley Seismic Zone and Imperial fault may be associated 
with the propagation of the Gulf of California rift system 
into the North American continent. 
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Joint Inversion of Shear Wave Travel Time Residuals and Geoid and 
Depth Anomalies for Long-Wavelength Variations in Upper Mantle 
Temperature and Composition Along the Mid-Atlantic Ridge 
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We report measurements of 55-5 differential travel time residuals for nearly 500 paths crossing the 
northern Mid-Atlantic Ridge. Differential travel times of such phases as SS and S with identical source 
and receiver have the advantage that residuals are likely to be dominated by contributions from the 
upper mantle near the surface bounce point of the reflected phase (55). Under this assumption, 
differential 55-5 travel time residuals are mapped at the 55 bounce points as a means of delineating 
lateral variations in mantle structure. After removing the signature of lithosphere age, we find evidence 
for variations in 55-5 residuals along the ridge at wavelengths of 1000-/000 km. These travel time 
anomalies correlate qualitatively with along-axis variations in bathymetry and geoid height. We 
formulate a joint inversion of travel time residual, geoid height, and bathymetry under the assumption 
that ail arise from variations in upper mantle temperature or bulk composition (parameterized in terms 
of Mg #). The inversion employs geoid and topography kernels which depend on the mantle viscosity 
structure. Inversion for thermal perturbations alone provides good fits to travel time and geoid data. 

The fit to topography, which is likely dominated by unmodeled crustal thickness variations, is not as 
good. The inversions for temperature favor the presence of a thin low- viscosity layer in the upper 
mantle and temperature perturbations concentrated at depths less than 300 km. Compositional 
variations alone are unable to match the travel time and geoid or bathymetry data simultaneously. A 
joint inversion for temperature and composition provides good fits to both geoid and travel time 
anomalies. Temperature variations are —50 K and compositional variations are ^0.5-3% Mg # for 
models with the temperature variations uniformly distributed over the uppermost 300 km and the 
compositional variations either distributed uniformly over the same interval or concentrated at 
shallower depths. The magnitudes of these variations are consistent with the chemistry and 
geothermometry of dredged peridotites along the Mid-Atlantic Ridge. 


Introduction 

Seismic velocity and density of upper mantle material are 
expected to be functions of temperature and composition. 
The delineation of long-wavelength variations in these phys- 
ical properties thus provides important constraints on mantle 
convection, crust-mantle differentiation, and mantle chemi- 
cal heterogeneity. In this study we determine lateral varia- 
tions in upper mantle temperature and composition along the 
Mid-Atlantic Ridge through the combined inversion of shear 
wave differential travel times, geoid height, and bathymetric 
depth anomalies. 

The advent of seismic tomography has led to a number of 
three-dimensional maps of lateral variations in seismic ve- 
locity in the upper mantle, and several such models of the 
North Atlantic region have been developed, both as parts of 
global studies [e.g., Woodhouse and Dziewonski , 1984; 
Nakanishi and Anderson , 1984; Tanimoto , 1990] and 
through regional investigations of long-period surface waves 
[e.g., Honda and Tanimoto , 1987; Mocquet et ai. f 1989; 
Mocquet and Romanowicz , 1990]. With surface wave meth- 
ods, each wave samples the average vertical variation in 
upper mantle structure along its path, but because of the long 
wavelengths involved, the inversion of phase or group 
velocity from many paths tends to smooth out lateral varia- 
tions. Body wave travel times can provide independent 
information about upper mantle heterogeneity at potentially 
shorter horizontal scales than surface waves can resolve, 
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and progress has been made in the determination of lateral 
heterogeneity in the North Atlantic through the use of both 
differential and absolute travel times of body waves [Kuo et 
a/., 1987; Grand, 1987, 1989]. 

The travel times used in this study are differential times of 
the body wave phase pair 55-5. Differential travel times of 
shear wave pairs are well suited to the study of upper mantle 
heterogeneity [Sipkin and Jordan , 1976, 1980; Stark and 
Forsyth , 1983; Butler, 1979; Kuo et al. , 1987; Woodward and 
Masters , 1991a] and have the advantage that source and 
receiver effects are approximately common to both phases 
and are thus largely eliminated by differencing. Under the 
assumption that the lower mantle is relatively homogeneous 
and that the portions of the wave paths in the upper mantle 
are steep, the differential travel time anomaly can be asso- 
ciated with upper mantle structure within a small volume 
centered beneath the surface bounce point of the reflected 
(55) phase. This technique is thus well suited to the inves- 
tigation of horizontal variations in structure, but the resolu- 
tion of variations with depth is poor. 

Oceanic bathymetry and geoid height data are sensitive to 
variations in mantle density at depth. Such variations can be 
either thermal or compositional in origin and, like seismic 
velocity, are presumably related to mantle convection and 
differentiation. Geoid (or gravity) and topography have 
become the most commonly used tools for mapping and 
constraining models of upper mantle convection [e.g., 
Anderson et a/., 1973; McKenzie and Bowin y 1976; McKen- 
zie , 1977; McKenzie et a/., 1980; Parsons and Daly , 1983; 
Buck and Parmentier, 1986; Craig and McKenzie , 1986]. In 
addition, measurement of the admittance (the spectral ratio 
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of geoid 10 topography) has been widely utilized to estimate 
the depth and mode of compensation of oceanic swells and 
plateaus [e.g., Warts et al., 1985; Cazenave et al . , 1988; 
Sandwell and MacKenzie , 1989; Sheehan and McNutt, 
1989]. Several workers [Dziewonski et al., 1977; Nakanishi 
and Anderson , 1984; Tanimoto and Anderson, 1984: Stark 
and Forsyth , 1983; Dz'iewonski, 1984; Kuo et al . , 1987] have 
noted correlations of geoid and travel time (or velocity 
structure) at a number of different wavelengths, although 
only a few [Hager et al., 1985; Hager and Clayton , 1989; 
Hager and Richards , 1989] have combined observational 
seismology with geoid anomalies in a quantitative and dy- 
namically consistent manner. 

Since very different convective flows can produce the 
same geoid and surface topography, the inversion of these 
data alone for the thermal or compositional source function 
is nonunique. Because this inverse problem is not well 
posed, most studies have concentrated on forward model- 
ling, i.e., varying the parameters of a simple model until a 
good fit to the data is achieved. With this approach, there is 
no guarantee that the set of parameters which give the best 
fit to the data is unique and that the correct solution has been 
isolated. Including seismic data provides additional con- 
straints which are sufficient to allow 7 us to formulate simple 
one-dimensional inversions. These inversions, of course, are 
still nonunique, but they are better constrained than those 
using only geoid or bathymetry data. 

In this paper we present the first formal inversion of geoid, 
depth, and 55-5 differentia! travel time anomaly data for 
lateral variations in upper mantle temperature and composi- 
tion along the Mid-Atlantic Ridge. Given a distribution of 
temperature or density perturbations in the upper mantle, 
the forward problem of calculating differential travel time, 
geoid, and depth anomalies is straightforward. This forward 
problem forms the basis for a joint linear inversion of these 
three types of observations under the assumption that all 
arise from parameterized long-wavelength variations in up- 
per mantle temperature or composition. Results of a set of 
inversions carried out under different assumptions regarding 
the depth extent of lateral heterogeneity and the mantle viscos- 
ity structure are compared with other constraints on variations 
in mantle temperature and degree of melt removal. 

Measurement of Differential Travel Times 

The seismic data used in this study consist of long-period 
5 and 55 phases obtained from the Global Digital Seismic 
Network (GDSN) [Peterson et al . , 1976; Peterson and Hutt , 
1982]; the Network of Autonomously Recording Seismo- 
graphs (NARS), a linear broadband array in western Europe 
[Nolet and Vlaar , 1982]; and several broadband stations 
from the global GEOSCOPE network [Romanowicz et al., 
1984, 1991]. We use only transversely polarized (SH) seis- 
mograms (rotated from N-S and E-W components) to avoid 
interference from the SKS phase and contamination from 
P-SV conversions at the base of the crust and other near- 
surface discontinuities. Recent work by Gee and Jordan 
[1989] suggests that travel times depend on the frequency 
band used in the analysis. In order to maintain a self- 
consistent data set for our study, additional processing is 
applied to data from the NARS and GEOSCOPE networks 
in order to mimic the instrument response of the longer- 
period GDSN stations. This processing allows us to measure 
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Fig. I. An example of the measurement of 55*5 differential 
travel time for the event of December 24, 1985 (10 km focal depth), 
at GDH (63° epicentral distance), (a) “Synthetic" 55 pulse gener- 
ated from 5. The 5 pulse is windowed and attenuated to account for 
the additional time that 55 travels in the mantle (r* = 3 s), and a tt!2 
phase shift is applied. ( b ) Windowed 55 wave pulse, (c) Cross 
correlation of the trace in Figure 1 b with that in Figure 1 g. The 
differential travel time residual is -5.0 s. 


travel times from a set of seismograms that all have essen- 
tially the same frequency response. Data from the NARS 
and GEOSCOPE stations are decimated (with a low-pass 
antialiasing filter) to a common sampling interval of 1 s. The 
data are further filtered using a noncausal three-point But- 
terworth filter [Rader and Gold , 1 967] with a frequency band 
pass of 0.01-0.20 Hz. This additional filtering greatly im- 
proves the signal-to-noise ratio of the 55 phase. 

A waveform cross-correlation method is utilized to deter- 
mine the differential travel time between the phases 5 and 
55 [ Butler , 1979; Stark and Forsyth, 1983; Kuo et al, 1987]. 
The procedure involves the construction of a “synthetic'’ 
55 pulse from 5 and the evaluation of the cross-correlation 
function between the real and synthetic windowed 55 
phases (Figure 1). The synthetic 55 pulse is created from 5 
in the following manner. The 5 pulse is windowed and 
attenuated (with attenuation parameter /* = 3 s) [Grand and 
Helmberger , 1984; Kuo et al., 1987] to account for the 
additional time 55 travels in the mantle, and then a rril 
phase shift (Hilbert transform) is applied to the attenuated 5 
pulse to simulate the frequency-dependent phase shift which 
the 55 wave undergoes at an internal caustic [Choy and 
Richards , 1975]. The differential time is obtained from the 
peak of the cross correlation between the synthetic 55 
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constructed from the 5 wave and the real 55. The residual 
^5-5 times are obtained by subtracting the observed differ- 
ential time from that predicted by the preliminary reference 
earth model (PREM) [Dziewonski and Anderson , 1981] and 
correcting for Earth ellipticity [Dziewonski and Gilbert , 1976] 
and 55 bounce point bathymetry. Our convention is that 
negative residuals are indicative of either early 55 or late 5. 

Constant window lengths of 120 s are used for both the 5 
and 55 phases. In general, the observed differential travel 
times vary by as much as 1 s depending on how 5 and 55 are 
windowed. Our modeling with synthetic seismograms indi- 
cates that emphasizing the onset of the 55 waveform can 
lead to bias for bounce points in areas of oceanic sediments. 
The effect of sediments at long periods is to produce precur- 
sory arrivals from reflections at the base of the sediments 
and late arrivals from waves which travel through the 
low-velocity sediments and are reflected at the crust-water 
interface. The net effect, after convolving the crustal re- 
sponse with the long-period GDSN instrument response, is 
that the time center of the 55 phase is effectively unchanged 
but the pulse is broadened both at the front and at the back. 
In our procedure the use of a constant window containing 
the entire 55 pulse should yield differential travel times that 
are little affected by the presence of sediments. 

Data 

The North Atlantic is an ideal area for conducting a 
differential travel time study in terms of the geographic 


distribution of available events and stations at suitable 
distances. The range in source-receiver separation was taken 
to be 55°— 86° to ensure separation of 5 and ScS at the longer 
distances and to avoid triplication in 55 at shorter distances. 
The 55 and 5 phases bottom from about 670 to 2300 km 
depth. We performed a search over all earthquakes in the 
Harvard centroid moment tensor (CMT) catalog for the 
years 1977-1987 [Dziewonski et al., 1981; Dziewonski and 
Woodhouse. 1983] and over all GDSN, NARS, and GEO- 
SCOPE digital seismic stations in order to find event-station 
pairs of the proper epicentral distance which provide 55 
bounce points in the North Atlantic region. Epicenters were 
obtained from the bulletin of the International Seismological 
Centre (ISC) for events occurring before 1987 and from the 
Preliminary Determination of Epicenters of the U.S. Na- 
tional Earthquake Information Service (NEIS) for events 
occurring in 1987. The final distribution of sources and 
stations used to measure 55-5 differential travel times is 
shown in Figure 2. The majority of data in this study comes 
from records of equatorial fracture zone earthquakes at 
North American and European stations, north and central 
Atlantic events at North American stations. Central Ameri- 
can events at European stations, and Mediterranean and 
European earthquakes at North American stations. 

This search yielded over 2000 event-station pairs with the 
proper epicentral separation. After winnowing the list be- 
cause of station inoperation, poor signal-to-noise ratio for 
the phases of interest, and interfering events, the final data 
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Fig. 3. (a) Map view of SS-S residuals relative to PREM [Dziewonski and Anderson, 1981], corrected for Earth 

ellipticity and seafloor bathymetry. Residuals are plotted at the 55 bounce point. The size of each symbol scales linearly 
with the magnitude of the residual. Lambert equal-area projection with pole of projection at 40°N, 40°W. Negative 
residuals indicate either early 55 or late S. Plate boundaries are from DeMets et al. [1990]. ( b ) Same as Figure 3 a but 
including data only for SS bounce points on lithosphere younger than 100 Ma. 


set consists of nearly 500 SS-S differential travel time 
residuals with bounce points in the North Atlantic (Figure 3). 
Uncertainties are determined for each measurement follow- 
ing the procedure outlined in Appendix A. A tabulation of all 
residuals, by station and source, is given by Sheehan [1991]. 

Results 

We interpret the variations in SS-S differential travel 
times in terms of lateral velocity variations within the crust 
and upper mantle beneath the surface reflection points of the 
SS wave path. Kuo et al. [1987] and Woodward and Masters 
[1991a] tested the validity of this assumption by plotting 
absolute S and SS residuals against SS-S residuals. They 
found that S and SS-S residuals are uncorrelated, while SS 
and SS-S residuals are strongly correlated, indicating that 
the assumption is justified. The validity of this assumption is 
further supported by the strong correlation of SS-S times 
with surface tectonic features in the vicinity of the SS 
bounce point. The residuals are further interpreted in terms 
of such upper mantle processes as lithospheric aging, flow- 
induced anisotropy, and along-axis heterogeneity in mantle 
structure. 

Lithospheric Aging 

Cooling and thickening of the lithosphere should yield a 
tendency toward an increase in seismic velocity with in- 


creasing lithospheric age. A linear regression experiment 
was performed to examine the correlation of the SS-S 
residuals with seafloor age. A gridded map of seafloor ages 
was constructed for the North Atlantic from the magnetic 
anomalies of Klitgord and Schouten [1986] and ages assigned 
according to Kent and Gradstein [1986] and Klitgord and 
Schouten [1986]. The isochrons of Sclater et al . [1981] were 
used in a few regions which were not covered by the 
Klitgord and Schouten [1986] data set. To obtain a represen- 
tative age value for the region spanning approximately one 
horizontal wavelength of the incident (55) wave, an average 
seafloor age was estimated for a 1° x 1° box centered on each 
55 bounce point. To reduce scatter, measurements whose 
bounce point depths differed by more than 2500 m from the 
depth predicted by the Parsons and Sclater [1977] plate 
cooling model were excluded from the final age regression. 
Although each 55 wave samples the upper mantle at a finite 
range of lithosphere ages, we expect that the different travel 
time anomalies contributed by the 55 path segments on the 
younger and older sides of the bounce point approximately 
cancel so that the age at the 55 bounce point is appropriate 
to the associated 55-5 residual. 

The 55-5 residuals for the North Atlantic are consistent 
with the expectation of an increase in seismic velocity with 
seafloor age. For bounce points between 0° and 60°N lati- 
tude, the coefficient derived by linear regression of residual 
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Fig. 4. 55-5 travel time residual versus square root of seafloor 

age for data from 0 to 60°N. Each plotted point represents the 
weighted mean of 14 adjacent data points. Weights are constructed 
from variances determined as discussed in Appendix A. Horizontal 
and vertical bars are ^standard errors of the means of the travel time 
residuals and (age) 1 ' 2 . Linear regression yields a slope of -0.68 := 
0.08 s Ma -,/2 fora 0-100 Ma age range (solid line) or -0.76 ± 0.09 
s Ma _1/ ‘ for a 0-80 Ma range (dashed line). 


v Th square root of age is -0.68 ± 0.08 s from 0 to 

\ j0 Ma, with a linear correlation coefficient of -0.85 (Figure 
4). Residuals north of 60°N, however, do not seem to be 
strongly correlated with lithospheric age. This may be due to 
the fact that this area is more tectonically complicated than 
"normal” oceanic lithosphere [e.g., White, 1988; Zehnder 
and Mutter , 1990], includes several ridge jumps, is in close 
proximity to continental regions and does not closely follow 
the age-depth relation of Parsons and Sclater [1977]. Com- 
pared with the residuals for 0-60°N, those from north of 
W*N are anomalously negative at young ages and anoma- 
lously positive at older ages. Although the age dependence 
of data from these high latitudes is anomalous, these data are 
no less reliable and will be retained in our analysis of 
heterogeneity. The slope of 55-5 residual versus square root 
of age for data from 0 to 60°N is smaller than that inferred 
from 5 delays of intraplate earthquakes in the Atlantic by 
Duschenes and Solomon [1977] (two-way S delay = — 1.2 s 
Ma“ l/2 ) and that reported by Kuo et al. [1987) (-1 s Ma~ l/2 ). 
It is larger, however, than the global average obtained by 
Woodward and Masters [1991a] (—0.5! s Ma~ l/2 ). Because 
me residual-age relation is not constant over the entire North 
Atlantic, some of these variations in slope may reflect real 
geographic differences. 

We may compare the variation of SS-S residual versus 
age with that due only to lithospheric cooling. For a litho- 
spheric structure given by the plate cooling model of Par- 


sons and Sclater [1977], we may convert temperature vari- 
ations to differences in shear velocity v s by adopting a value 
for dv s ldT , which we take to be uniform and equal to -0.6 
m/s K" 1 [McNutt and Judge , 1990]. For a horizontal slow- 
ness typical of the teleseismic 5 and SS waves of this study 
(0.1375 s km” 1 ), the slope of the line best fitting the 55-5 
travel time delay versus age given by the plate cooling model 
over 0-100 Ma is then —0.64 ± 0.07 s Ma~ ,/2 , a result 
indistinguishable from the observed slope. This agreement 
indicates that the dependence of travel time residual on plate 
age can be explained entirely by lithospheric cooling. 

The trend of the travel time residual versus lithospheric 
age relation changes at about 100 Ma. After 100 Ma, the 
residuals appear to flatten out (Figure 4), in the same sense 
as the plate cooling model of Parsons and Sclater [1977]. 
Such a pattern may reflect the unmodeled effect of increased 
sediment or crustal thickness or, as suggested by Parsons 
and Sclater [1977] may be partially the result of secondary 
convection which supplies heat to the base of the plate at 
older ages. To avoid possible biases associated with any of 
these effects, we shall restrict our analysis to data with 
bounce points on seafloor younger than 100 Ma. To look for 
other systematic variations in the residuals, we correct for 
age by removing the linear relation shown by the solid line in 
Figure 4. This correction is effectively a normalization of 
residuals to 22-Ma-oId lithosphere (the zero crossing of the 
regression line). Since the age correction is linear, normal- 
ization of the residuals to a different age would simply result 
in a constant value being added to or subtracted from all 
residuals. 


Anisotropy 

Another systematic velocity variation that has been sug- 
gested as a possible contributor to residual 55-5 travel times 
is azimuthal anisotropy. Kuo et al. [1987] examined this 
phenomenon in detail and concluded that alignment of 
olivine crystals in the asthenosphere created a significant 
pattern of azimuthal anisotropy in 55-5 residuals measured 
in the Atlantic region. We have also searched for evidence of 
azimuthal anisotropy with our data set. 

Backus [1965] and Crampin [1977] demonstrated, from the 
general form of body wave anisotropy in a weakly anisotro- 
pic medium, that the linear form of the azimuthal variation of 
velocity is given by 

V 2 = A 0 + A| cos 28 + A 2 sin 26 4* Ay cos 4 8 

4- A 4 sin 4<9 (1) 

where V is the body wave velocity, the A n are linear functions 
of the elastic moduli, and 8 is an azimuth, defined for our 
problem by the angle between the great circle path and the 
direction to geographic north measured at the 55 bounce point. 
Equation (1) was further simplified by Kuo et al. [1987] and 
parameterized in terms of travel time residuals: 

R = R q + /?, cos 28 + R z sin 28 

+ R 3 cos 4 8 4* R a sin 4 8 (2) 

where R is the travel time residual and the R n are constants. 
By fitting a function of this form to our age-corrected 
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Fig. 5. Age-corrected 55-5 residual (see text) versus azimuth 0 . 
Each plotted point represents the weighted mean of 10 adjacent data 
points. The solid curve shows the best fitting 40 variation derived 
from these data. The dashed curve shows the preferred model of 
Kuo ei al. [1987], which corresponds to an alignment of the a axis of 
olivine in the approximate direction N I3°W. 


measurements we can determine if our data are consistent 
with the presence of anisotropy. 

We have conducted several tests of azimuthal anisotropy 
with our travel time data. We performed least squares 
inversions to determine 2 0 and 40 patterns which provide 
best fits to the age-corrected 55-5 residuals. Data from the 
region north of 60°N were included in the regressions dis- 
cussed here, but a second set of regressions performed 
without these data yielded similar results. The anisotropy 
indicated by our regression experiments differs significantly 
from the preferred model of Kuo et al. [1987] both in 
magnitude and in phase (Figure 5). Our results indicate that 
for the 26 model the slow direction for 55-5 is N4°W and the 
peak-to-peak magnitude of the effect is less than 1 s; for the 
46 model the slow directions are N32°W and N58°E and the 
magnitude is 2.5 s; for the joint 20 and 40 model the slow 
direction is N32°W and the magnitude is just under 3 s. Kuo 
et al. [1987] obtained a peak-to-peak variation with azimuth 
of 5-7 s and a slow direction at N13°W. The slowest 
residuals in the Kuo ex al. [1987] study were from north- 
south paths, i.e., nearly along the ridge, and the fastest 
residuals were from northeast-southwest trending paths with 
bounce points north of the Azores-Gibraltar plate boundary 
(an area noted to be anomalously fast in their study), so their 
reported anisotropy may have been at least partly the result 
of unmodelled upper mantle heterogeneity. Our inversion for 
a 20 pattern of anisotropy provided a variance reduction of 
only 2%, compared with 20% for a 40 pattern and 22% for a 
combined 20 and 40 pattern. On the basis of these values of 
variance reduction and the number of free parameters in- 
volved, our results suggest that there is no single coherent 
pattern of upper mantle anisotropy in the North Atlantic. 
Inversions were performed on geographic subsets of the data 
to test whether there may be several distinct patterns of 


anisotropy acting in the North Atlantic region, but the 
results were inconclusive because of a poor sampling of 
azimuths in the smaller subregions. The latest anisotropic 
upper mantle models obtained from surface wave tomogra- 
phy [Montagner and Tanimoio , 1990] show a complex 
pattern of anisotropy in the North Atlantic region. Any 
azimuthal anisotropy in the asthenosphere induced by plate 
motions in the North Atlantic may be heterogeneous be- 
cause the three plates in the region are slow moving and the 
return flow is not closely related to plate divergence [Hager 
and O'Connell , 1979, 1981; Parmentier and Oliver , 1979]. 

Spatial Patterns of Age-Corrected Residuals 

After removal of the dependence on seafloor age, a plot of 
55-5 travel time residuals at the 55 surface reflection point 
(Figure 6) shows several interesting features. Perhaps the 
most striking is that residuals in the western Atlantic north of 
about 35°N are on average nearly 4 s more negative than 
those to the south. This feature is also noticeable in Figure 3 
but is more obvious after age dependence is removed. A 
similar change at approximately this latitude was noted for 
55-5 residuals with 55 bounce positions on the eastern side 
of the Mid-Atlantic Ridge by Kuo et al. [1987] and was 
attributed to a change in upper mantle structure across the 
Azores-Gibraltar plate boundary. The signal that we observe 
is predominantly from data with bounce points on the 
western side of the ridge. A map view' of the azimuthal 
distribution is shown in Figure 7 and serves as an aid to 
assess qualitatively the geometry of wave paths to the south 
and north of 35°N. We examined the possibility that this 
signal may be from the Caribbean anomaly, a region of 
anomalously high velocity in the mantle between 600 and 
1400 km depth beneath the Caribbean originally reported by 
Jordan and Lynn [1974] and further confirmed by Grand 
[1987], If the first leg of the 55 wave paths from South 
America to western Europe were to bottom in the high- 
velocity Caribbean region, the result would be early 55-5 
residuals. This would produce a feature of opposite sign 
from that observed, so we discount it as an influence here. 
Another possible explanation for the long-wavelength signal 
could be azimuthal anisotropy, but the examination above of 
possible patterns of azimuthal anisotropy does not support 
this suggestion. 

Another distinctive feature of the residuals in Figure 6 is a 
row of negative values which trends northwest to southeast 
along the trace of the New England Seamounts and across 
the ridge to the vicinity of the Great Meteor Seamount. This 
feature in the residual map comes from event-station pairs at 
a number of different azimuths and distances, so cannot be 
attributed to a source or receiver effect. We do not observe 
distinctive anomalies in the vicinity of the Bermuda, Azores, 
or Canary Islands hotspots. The data density is poor for the 
Bermuda region, however, and any signal associated with 
the Canary Islands may be obscured by the ocean-continent 
transition. Recently active hotspot islands might be ex- 
pected to display strong positive (late) residuals, such as 
Stewart and Keen [1978] observed for PP-P residuals at the 
Fogo Seamounts. In contrast, Woodward and Masters 
[1991a] found mostly negative (early) SS-S residuals in the 
vicinity of the Hawaiian hotspot, and Jordan [1979] and 
Sipkin and Jordan [1980] have suggested that the net effect 
of hotspots may be to produce early arrivals because of the 
presence of high velocities in a depleted mantle residuum. 
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Fig. 6. (a) Map view of age-corrected SS-S residuals. ( b ) Same as Figure 6 a but including data only with 55 bounce 

points on lithosphere younger than 100 Ma. 


There is a systematic variation of 55-5 residual with 
latitude, i.e., effectively along the direction of the Mid- 
Atlantic Ridge axis. Age-corrected 55-5 residuals with 55 
bounce points on lithosphere younger than 100 Ma are 
shown versus latitude in Figure 8. The along-axis variations 
show a variety of scales, notably at wavelengths of about 
1000-2000 km in the region from 15° to 35°N, and at about 
6000 km wavelength from late (positive residuals) in the 
south (20°-35°) to early (negative residuals) farther north 
(45°-55°N). The largest of these variations are robust with 
respect to selective removal of portions of the data. The 
Iceland region appears as a local maximum (positive 55-5 
delay) on the profile, but the Azores hotspot does not have a 
distinct seismic signal. 

Joint Inversion of Travel Time Residuals 
and Geoid and Depth Anomalies 

Long-wavelength variations in shear wave velocity of the 
sort depicted in Figure 8 presumably are a consequence of 
some combination of variations in temperature and compo- 
sition of the upper mantle. Such lateral variations should 
also have signatures in other physical quantities measurable 
at these wavelengths, notably gravity (or geoid height) and 
topography (or residual bathymetry), because of the depen- 
dence of these quantities on bulk density. Travel time 
residuals, geoid anomalies, and residual depth anomalies are 
independent quantities dependent in different ways on tem- 
perature, bulk composition, and their variation with depth. 
We therefore seek a quantitative procedure for treating 


travel time residuals jointly with geoid and bathymetry data 
and in particular for a combined inversion of all three 
quantities for horizontal variations in upper mantle temper- 
ature and composition. 

To ensure complementarity of data sets, bathymetry and 
geoid height values are obtained at each 55 bounce point, 
and both are corrected for subsidence with seafloor age by 
means of the plate cooling model [Parsons and Sclater , 1977; 
Parsons and Richter , 1980]. In this manner we effectively 
normalize all observations to zero age. Bathymetric data are 
obtained from the corrected Digital Bathymetric Data Base 
(5 arc min grid) [U.S. Naval Oceanographic Office , 1985]. 
Geoid data are taken from a combined set of Seasat and 
GEOS 3 altimeter data [Marsh et al 1986]. Data north of 
70°N were not included in the Marsh et al. [1986] data set 
due to the high probability of being over sea ice, so our 
analysis below is confined to latitudes less than 70°N. We 
find that the correlation of 55-5 residuals with the low-order 
geoid is negative but that at high order the correlation is 
positive (Figure 9). This relationship may indicate a depth 
dependence of contributions to geoid and travel time (e.g., 
the long-wavelength signal may be a lower mantle effect). 
Low-degree harmonics are likely linked to deep-seated den- 
sity heterogeneities and subducting slabs [Hager, 1984; 
Hager et al ., 1985]. Since we are interested in upper mantle 
processes, we filter out the long-wavelength component of 
the geoid by subtracting a reference field [Lerch et al., 1979] 
expanded in spherical harmonics to degree and order 7 and 
tapered to degree and order 11. To provide a comparable 
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Fig. 7. Map view of the distribution of sampling azimuths. Lines indicate the wave path azimuth at the 55 bounce 

point. Mercator projection. 


bathymetric data set, bathymetry is high-pass-filtered (cor- 
ner at 4000 km, cutoff at 6000 km) to remove long- 
wavelength trends. Aiong-axis profiles are constructed from 
the age-corrected and filtered geoid and bathymetry data. 
These profiles, of course, are only approximations to true 


along-axis variations, but they share many features with 
similarly filtered geoid and bathymetry profiles obtained 
from observations along the ridge axis rather than at the SS 
bounce point locations. 

Profiles of age-corrected travel time residuals, geoid, and 
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Fig 8 Age-corrected 55-5 residual versus latitude along the Mid-Atlantic Ridge frorn 10° to 90°N . The residuals 
shown are moving averages (such that each point is used twice) of 12 adjacent data points. Error bars indicate standard 
errors of the mean residual and position. Bounce points on lithosphere of age 0-100 Ma are used The approximate 
locations of several fracture zones (Fifteen-Twenty. Kane. Atlantis. Oceanographer, and Charlie-G.bbs. denoted by 
abbreviations) and of the Iceland and Azores hotspots are indicated. 
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Fig. 9. Linear correlation, by highest harmonic degree removed 
from the geoid, of observed 55-5 residual with geoid height 
measured at the corresponding 55 bounce point. Both travel lime 
and geoid residuals are age-corrected. First, the uncorrccted geoid 
data [Marsh et al., 1986] are correlated with 55-5 residuals and a 
slope and correlation coefficient determined. Then a geoid reference 
field [Lerch et a/., 1979] up to degree and order 2 (with taper to 
degree and order 6) is calculated and removed from the geoid data, 
the slope and correlation coefficient with 55-5 residuals calculated, 
and so on for higher harmonic degrees t removed from the geoid 
data, with appropriate tapers (up to t + 4). (a) Linear correlation 
coefficient between geoid and 55-5 residuals versus highest har- 
monic degree and order removed from the geoid. (b) Slope of the 
correlation between geoid and 55-5 residual data, as a function of 
highest harmonic degree and order removed from the geoid. Extra 
points at degree and order 10 are obtained by using different tapers 
(no taper, taper to i - 14, and taper to i - 15). 


bathymetry are compared in Figure 10. All profiles have 
been smoothed by applying a running average. While quali- 
tative correlations among profiles are apparent, we seek to 
quantify possible models of temperature and compositional 
variations that can match these observations. Oceanic 
bathymetry and geoid height are both sensitive to variations 
in mantle density at depth. Such variations can be either 
thermal or compositional in origin and, like seismic velocity, 


SS-S RESIDUAL 




Fig. 10. Comparative plots of age-corrected (<a) 55-5 residual, 
(b) geoid, and (c) bathymetry, along the Mid- Atlantic Ridge, 10°- 
65°N. Bathymetry and geoid have been high-pass- filtered (see text). 
All of the data points shown are moving averages of 10 adjacent data 
points. Bounce points from Lithosphere of age 0-100 Ma are used, 
except that data from the Labrador Sea region are omitted. The solid 
line represents a filtered version of the depicted observations, contain- 
ing only the wavelengths used in the inversions (1400-7100 km). 


are presumably related to mantle convection and differenti- 
ation. For a given density change, the seismic signature of 
thermal and compositional heterogeneity are of opposite 
sign, so travel time residuals constitute key information for 
distinguishing between mechanisms of heterogeneity. 

Inversion for Thermal Structure 

We seek to formulate an inversion for the distribution of 
temperature anomalies T( x, z) (where x is along-axis and z is 
depth) that can produce the along-axis geoid, bathymetry, and 
travel time anomalies shown in Figure 10. Topography and 
geoid kernels were calculated for prescribed models of viscos- 
ity for an incompressible, self-gravitating, Newtonian mantle 
with free slip at the surface and the core-mantle boundary. The 
convecting region is assumed to be overlain by a high-viscosity 
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layer 40 km thick. We performed calculations both for a mantle 
of constant viscosity and for a mantle with a shallow low- 
viscosity layer. Topography and geoid anomalies depend on 
the viscosity structure, but the predicted travel times do not. 
Kernels were calculated using a method similar to that of 
Richards and Hager [1984] except that the solution was di- 
rectly integrated across the layers instead of being obtained via 
propagator matrices [McNutt and Judge, 1990]. 

The inversion is best conducted in the horizontal wave- 
number domain. The thermal anomalies A7X*, z) at depth 
are related to the predicted dynamic topography h(k) for 
wavenumber k via an integral of the form 


AM*) = — P“‘ H(k , z)ATT*, z) dz (3) 

P o “ P». J, 

[Parsons and Daly , 1983], where a is the volumetric coeffi- 
cient of thermal expansion, p 0 and p w are the densities of the 
mantle and of water at standard temperature and pressure, 
and £ min an ^ Zmax are the upper and lower boundaries of the 
layer in which temperatures are allowed to vary. Table 1 
contains a summary of the constants adopted here. The 
depth- and wavenumber-dependent topography kernel H{k, 
z) is calculated from the equations of continuity and motion 
given a set of boundary conditions, a viscosity model, and a 
constitutive relation between stress and strain [Parsons and 
Daly , 1983]. Similarly, the kernel G(* , z) for the geoid 
relates the thermal anomalies to the geoid N(k) via 

2 ffTpQa f Zm»x 

A N{k) = 7 G(k, z)A7(*, z) dz (4) 

Qk J . 

[Parsons and Daly , 1983], where T is the gravitational 
constant and g is the surface gravitational acceleration. 

Sample geoid and topography kernels calculated for dif- 
ferent wavenumbers and viscosity structures are shown in 
Figures 1 1 and 12. Cartesian kernels are used throughout this 
study because of their computational efficiency and straight- 
forward application to Fourier transform techniques. We 
have compared extrema of the upper mantle portions of the 
geoid and topography kernels for a layered Cartesian Earth 
and a spherical Earth for a number of wavelengths and 
different viscosity structures (Figure 12), and we note good 
agreement even at very long wavelengths (spherical har- 
monic order i = 6). This agreement suggests that the results 
presented here should be applicable to the spherical Earth 
without introducing unreasonably large errors. 

Temperature perturbations at depth can be converted to 
seismic velocity perturbations by assuming a value for the 
partial derivative of shear wave velocity with respect to 
temperature, BvjBT. The resulting two-way travel time 
perturbation is given by 


a m = 


= 2 — f Zm 

ST 

J *-mi 


ATT*, z) dz 
W,(Z) 2 [1 ~p 2 v,U) 2 } m 


(5) 


where v s (z) is from the reference shear velocity model 
[Dziewonski and Anderson, 1981] and p is the ray parameter, 
generally taken to be a representative value for the range of 
epicentral distances considered here. We use a value of -0.6 
ms -1 K" 1 for BvJBT. This value is more negative than the 
values of Anderson et al. [1968] and Kumazawa and Ander- 
son [1969] at standard temperature and pressure but is 


similar to the value of -0.62 m s" 1 K" 1 determined by 
McNutt and Judge [1990] by a least squares fit of Love wave 
phase velocities to predicted temperature of the lithosphere. 
Such a value is consistent with the change in P wave velocity 
with temperature, Bv p lBT = -0.5 m s“ l K“\ found from 
modeling wave propagation along subducting slabs [Creager 
and Jordan , 1986; Fischer et aL , 1988] if we assume that 
BvjBT = 1.1 dv p lBT [Woodhouse and Dziewonski , 1984]. 
Partial melt would increase the value of BvjBT [Sleep, 1974; 
Sato and Sacks , 1989], but simultaneous analysis of both 
shear and compressional differential travel times by Wood- 
ward and Masters [1991a] indicates that significant partial 
melting is not required to explain the differential travel time 
residuals in the North Atlantic region. 

The forward problem consists of calculating geoid, topog- 
raphy, and travel time residual profiles given a starting 
two-dimensional temperature structure T(x, z). The inverse 
problem consists of finding a temperature structure that 
predicts (via equations (3M5)) geoid, topography, and travel 
time profiles which best fit those observed. The familiar 
matrix equation d = A m is formed from discrete versions of 
equations (3)-(5). The data vector d consists of the topogra- 
phy, geoid, and travel time residuals, the model vector m 
contains the temperature variations for which we are solv- 
ing, and the matrix A contains the coefficients and kernels 
which relate the data to the model. As a check on ou; 
procedure, we constructed a forward problem for geoid and 
topography and found good agreement with the modelling 
results of McKenzie et al. [1980]. 

The observed profiles of bathymetry, geoid, and travel 
time are interpolated to a constant spacing, demeaned, 
tapered at both ends with a 10% sine squared taper, and 
Fourier transformed (Figure 10). Since the profiles for which 
these operations were performed extend from 10° to 72°N, 
the first and last 10% of the profile (10°-16°N and 66°-72°Nl 
will be affected by the taper. The 3 k x 1 data vector d is the . 
constructed, using the complex (to retain both amplitude and 
phase) bathymetry, geoid, and travel time data sampled at n 
discrete wavenumbers: 

d = [>/!(*,), ••• . A *(*„), 

A N(A-,), , A N(k„), Ar(*,), , M(k„)] T (6) 

where T denotes transpose and n in our problem is equal to 
5, representing the first five coefficients of the Fourier ser *s 
expansion (wavelengths of 7104, 3552, 2368, 1776, and 1420 
km). For the case where temperature perturbations are 
constrained to be in a single layer, the n x 1 model vector m 
is given by 


m = [ATT*,), ••• , A7-(*„)] r (7) 

For the more general case of a multilayer system, the nj x 1 
model vector m is given by 

m = [ATX*,, Z\), • • - , ATX* , , Z;),AT(* 2 , Zj), • ’ ’ - 

ATT * 2 . zj). • • • , ATT**, Z|), • • • , ATT*,,, ij)) T (8) 

where z< is the layer index and j is the total number of layers. 
In this paper we perform inversions for single-layer models 
only. The “layers” of temperature variations are independent 
of the “layering” system of lid, low-viscosity zone, and mantle 
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TABLE l. Adopted Constants 

Variable Description 

volumetric coefficient of thermal expansion 
average mantle density 
density of seawater 
gravitational constant 
surface gravitational acceleration 
thermal coefficient of shear velocity 
variation of shear velocity with Mg # 
variation of density with Mg # 
average SS ray parameter at 70° 

*Stacey [1977] and Duffy and Anderson [19891. 
t Akimoto [1972]. 


a 

Po 

Pw 

r 

9 

dvJdT 

dVj/dMg 

dp/d Mg 


Value 


2.5 x 10' 5 K - '* 

3300 kg m~ J 
1000 kg m~ 3 

6.67 x 10"" N m : kg' 2 

9.8 m s' 2 

-6.0 x 10~* km s' 1 K 1 

1.8 x 10' 2 km s' 1 Mg #"'t 
-12 kg m~ 3 Mg #~'t 
0.1375 s km -1 


(a) CONSTANT 


GEOID KERNEL 



VISCOSITY 



(b) LOW VISCOSITY LAYER 


GEOID KERNEL 



TOPOGRAPHY KERNEL 



Fig. 11. Upper mantle portion of the kernels for geoid and topography at two wavelengths A - 2 jrik for two 
viscosity models. The convecting region in both models is overlain by a h.gh-viscos.ty layer 40 km thick wuhviscoity 
10 4 that of the underlying mantle («) High-viscosity lid is underlain by a mantle ofumfonn v.scosity and other physic^ 
parameters, (b) High-viscosity lid is underlain by a zone extending to a depth of 200 km having a viscosity equal to 0.01 
that of the underlying mantle. 






19,992 


Sheehan and Solomon: Joint Inversion of Travel Time, Geoid, and Depth 


which we use for the calculation of kernels, although major 
changes in viscosity would tend to segment A T as well. The 
temperature layering simply refers to that region bounded by 
z m ^ and z max in the integrals of equations (1H 3 )* 

The 3 n x n matrix A contains the coefficients and kernels 
that relate the temperature perturbations to the observa- 
tions, which for the single-layer case is given by 


The resulting model vector m is inverse Fourier transformed 
back to the spatial domain to produce an along-axts temper- 
ature profile. The solution m is also substituted into equa- 
tions (3H5) to compare predicted geoid, bathymetry, and 
travel time residuals with those observed. 

To assess how well the various data are being fit, reduced 
X 2 values are calculated in the spatial domain from 




The matrix A contains both bathymetry and topography 
kernels and is thus viscosity dependent; that is, a viscosity 
structure must be assumed. We solve the equation d = A m 
by least squares 

m = (AR^r'AR^’d (10) 

where A is the complex conjugate transpose of the matrix A. 
Construction of the data covariance matrix Rjj is discussed 
in Appendix B. Equation (10) is solved for the solution 
vector m, and variance reduction is calculated via 

(d^lR^d-Am) 
variance reduction = 1 -r ryr 

d Rdd “ 


1 

2 t ~ y ?) 2 


where v = N - p is the number of degrees of freedom, N is 
the number of data points, p is the number of parameters 
being fit, a} is the estimated variance of the data (see 
Appendix B), and y p and y Q are the predicted and observed 
values of the smoothed profiles, respectively. If the model is 
a good fit to the observations, then xl should be approxi- 
mately unity. A xl greater than 1 indicates a poorer fit. A 
value of xl ! ess than * does not necessarily indicate a 
physically meaningful improvement in fit, however; it can 


(a) CONSTANT VISCOSITY (A = 4000 km) (c) CONSTANT VISCOSITY (A = 6667 km) 

GEOID KERNEL TOPOGRAPHY KERNEL GEOID KERNEL TOPOGRAPHY KERNEL 

-Q.4 OO Q.4 -0,4 0.0 0 4 0.6 -0,4 0.0 0.4 -0,4 0.0 04 0 8 
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TABLE 2. Inversion Models and Measures of Fit 


Layer 

Thickness, 

km 

Viscosity 

Structure 41 

ST 

Range, K 



Variance Reduction, % 



Reduced x 2 

Values 


aMg # 

Range Total 

Bathymetry Geoid 

SS-S 

Total 

Bathymetry 

Geoid 

SS-S 

0-150 

cvm 

180 


Temperaiure Variations Only 
53 25 

79 

58 

1.82 

2.90 

0.72 

1.84 

0-150 

lvz 

230 


57 

27 

79 

66 

1.68 

2.83 

0.72 

1.49 

0-300 

cvm 

60 


47 

21 

85 

41 

2.05 

3.09 

0.49 

2.57 

0-300 

lvz 

110 


57 

24 

85 

65 

1.65 

2.94 

0.49 

1.51 

0-650 

cvm 

20 


41 

14 

91 

25 

2.30 

3.33 

0.29 

3.27 

0-650 

lvz 

33 


49 

17 

83 

51 

1.96 

3.21 

0.55 

2.11 

0-150 

cvm 


1.1 

Compositional Variations Only 
33 46 

74 

-9 

2.60 

2.12 

0.89 

4.80 

0-150 

lvz 


2.4 

44 

75 

76 

-9 

2.19 

0.98 

0.81 

4.76 

0-300 

cvm 


0.4 

33 

29 

87 

-6 

2.62 

2.78 

0.43 

4.65 

0-300 

lvz 


1.3 

43 

73 

73 

-7 

2.22 

1.03 

0.92 

4.71 

0-650 

cvm 


0.1 

32 

19 

93 

-3 

2.63 

3.16 

0.22 

4.52 

0-650 

lvz 


0.5 

49 

65 

86 

+ 5 

1.99 

1.33 

0.46 

4.17 

0-150 

cvm 

210 

Thermal and Compositional Variations in Same Layer 
1.5 75 44 78 100 

1.01 

2.24 

0.78 

0.00 

0-150 

lvz 

235 

2.1 

86 

75 

80 

100 

0.56 

1.01 

0.69 

0.00 

0-300 

cvm 

110 

0.7 

73 

28 

89 

100 

1.09 

2.88 

0.39 

0.00 

0-300 

lvz 

125 

1.1 

84 

74 

76 

100 

0.62 

1.03 

0.83 

0.00 

0-650 

cvm 

55 

0.4 

71 

18 

94 

100 

1.15 

3.27 

0.20 

0.00 

0-650 

lvz 

60 

0.8 

85 

66 

88 

100 

0.58 

1.33 

0.42 

0.00 

0-150 

cvm 

Thermal Inversion in Lavers as Noted, Compositional Variations 0—50 km Only 

210 5.5 84 83 77 91 0.63 0.66 

0.81 

0.42 

0-150 

lvz 

240 

4.5 

85 

85 

70 

96 

0.60 

0.60 

1.04 

0.16 

0-300 

cvm 

80 

5.9 

80 

84 

86 

72 

0.80 

0.64 

0.50 

1.27 

0-300 

lvz 

120 

4.7 

86 

90 

77 

89 

0.56 

0.39 

0.80 

0.48 

0-650 

cvm 

25 

6.0 

73 

85 

92 

47 

1.07 

0.59 

0.27 

2.35 

0-650 

lvz 

35 

4.6 

75 

82 

73 

71 

0.98 

0.72 

0.95 

1.29 


*cvm, constant viscosity mantle; lvz, mantle with low-viscosity zone. 


simply be a consequence of overfitting data more closely 
than the estimated variance, xl values were calculated for all 
three data types simultaneously as well as for each individual 
data type. In the estimation of overall xl ♦ N is equal to the 
sum of the number of data in the three different profiles 
(residual geoid and bathymetry and differential travel time). 
In the estimation of xl values for individual profiles, N is the 
number of data points and p is equal to one third the number 
of parameters used in the inversion, since the xl value is 
being calculated using only one third the data going into the 
inversion. 

Six inversion experiments for temperature structure were 
performed (Table 2). Inversions were carried out for two 
different viscosity models and for three different thicknesses 
of the layer in which lateral temperature variations were 
assumed to occur. Because topography and geoid anomalies 
depend only on the ratios of viscosity in different layers 
[Richards and Hager , 1984; Robinson et ai , 1987; Hong et 
al., 1990], we set the dimensionless viscosity of the layer 
representing the bulk of the mantle to unity. In one viscosity 
model, termed the “constant viscosity mantle,” a 40-km- 
thick high-viscosity lid overlies a unit viscosity mantle. We 
set the viscosity in the lid to 10 4 , which effectively mimics 
rigid behavior. In a second model, a 160-km-thick low- 
viscosity zone is present beneath a 40-km-thick lid; the 
viscosity in the low-viscosity zone is a factor of 100 less than 
in the underlying mantle. The thickness of the layer of 
temperature perturbations was taken variously to extend 
from 0 to 150 km depth, 0 to 300 km depth, and 0 to 650 km 
depth. The matrix A is different for each of these cases, as it 


involves viscosity-dependent geoid and topography kernels 
and also a summation over depth. Since the geoid and 
bathymetry kernels do not include any phase information 
(except that a sign change produces a 180° phase shift), large 
differences in phase between the observed profiles of geoid 
and bathymetry would indicate that these components can- 
not simultaneously be well fit by our models. Spectral 
coefficients for both the filtered profiles and models are given 
by Sheehan [1991] and serve to indicate how well the various 
spectral components of the data are fit. 

Inversion results for the constant-viscosity-mantle cases 
are shown in Figure 13. Predicted profiles for each inversion 
solution were calculated from equation (5). For these solu- 
tions, the long-wavelength fit to the geoid is better than ai 
short wavelengths. The fit to bathymetry is poor. The 
predicted magnitude of the SS-S residuals range from a 
factor of 5 too small for the 650-km-thick layer to a factor of 
about 1.5 too small for the 150-km-thick layer. Increasing the 
temperature variations to improve the fit to the SS-S resid- 
uals leads to predicted geoid variations that are too large. 
The highest total variance reduction and best fit for the 
constant-viscosity cases come when lateral variations are 
constrained to shallow (0—1 50 km) depth. The variance 
reduction is 25% for bathymetry, 79% for geoid, and 58% fc; 
travel times, and the reduced x 2 values are 2.9 for bathym- 
etry, 0.7 for geoid, and 1.8 for travel times (Table 2). The 
total variance reduction is 53%, and the total xl value is 1 8. 
The variation in temperature is 180 K for the 150-km-thick 
layer and only 60 K for the 300-km-thick layer. 

Figure 14 shows inversion results for the models with a 
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Fig. 13. Results of combined inversion of geoid, bathymetry, 
and 55-5 travel time residuals for upper mantle temperature varia- 
tions. The viscosity structure is taken to consist of a 40-km-thick 
high-viscosity lid overlying a constant-viscosity half-space, (a) 
Three solutions for aiong-axis temperature variations: Dotted line 
indicates temperature perturbations constrained to be uniform over 
0-150 km depth. Long-dashed line indicates temperature perturba- 
tions constrained to be uniform over 0-300 km depth. Short-dashed 
line indicates temperature perturbations constrained to be uniform 
over 0-650 km depth. ( b ) Observed (solid line) and predicted 
aiong-axis profiles of 55-5 travel time residual. The “observed" profile 
is the filtered profile of Figure 10 and contains only the wavelengths 
used in the inversion (1400-7100 km). Predicted profiles were calcu- 
lated from equation (5). Line types correspond to those of the temper- 
ature models, (c) Observed and predicted aiong-axis geoid profiles. 
Same treatment as in Figure 13 b. (d) Observed and predicted aiong- 
axis bathymetry profiles. Same treatment as in Figure 1 3 b. 


thin low-viscosity zone. A good fit to both geotd and travel 
time is found, although the alignment in phase of predicted 
and observed geoid is not as good as for the constant- 
viscosity case. The fit to bathymetry is again poor. The total 
variance reduction for the 150-km-thick and 300-km-thick 
layers are both 57% (*; values are 1.68 and 1.65 respec- 
tively, so the total fit is slightly better for the 300-km-thick 
layer), although the shallow model provides slightly higher 
variance reduction for bathymetry (27% for 0-150 km deep 
layer 24% for 0-300 km deep layer: xl values of 2.8 and -9, 
respectively) and the 300-km-thick-layer model provides 
higher variance reduction for geoid (79% for 0-150 km deep 
layer. 85% for 0-300 km deep layer: xi- values of 0.7 and 0.5, 
respectively). The variation in temperature for the 150-km- 
thick layer is 230 K and in the 300-km-thick layer is 1 10 K. 

We have explored the hypothesis that the lack of correla- 
tion of predicted and observed topography is an indication 
that the source of variations in the geoid and travel time 
anomalies is deep. To test this hypothesis, we performed 
inversions with temperature variations restricted to deeper 
layers and found that fits to topography were still poor. It is 
possible that the bathymetric signal is dominated by crustal 
thickness variations which are not included in our calcula- 
tion of dynamic topography. An assessment of such thick- 
ness variations is discussed further in Appendix B. 


Inversion for Compositional Variations 

A possible alternative to aiong-axis variation in mantle 
temperature is lateral variation in bulk mantle composition, 
due perhaps to a variable extent of melt extraction or 
different degrees of mixing of compositionally distinct vol- 
umes of mantle material. The dynamical effects of composi- 
tionally induced density variations can be large [ O'Hara , 
1975; Bovd and McCallister. 1976: Oxburgh and Parmentier, 
1977: Sotin and Parmentier. 1989]. The fraction of mantle 
potentially extractable as basaltic melt is thought to be 
15-25% [e.g ..Green and Liebermann. 1976]. Thus, for every 
volume of basalt removed from the mantle, a volume of 
residuum several times larger is left behind. The effect _°f 
basalt depletion is to increase the molar ratio Mg/(Mg - re) 
(or Mg #) in the residuum, which reduces the density an 
increases the seismic velocities [e.g., Liebermann , 1970; 
Akimoto, 1972]. For example, subtraction of 20 mol % 
olivine basalt from pyrolite can decrease the density of the 
residuum by nearly 2%, equivalent to a thermal perturbation 
of nearly 500 K [ Jordan , 1979]. Thus compositional changes 
need only be slight to produce effects of the order of 100 K, 
comparable to values obtained from the inversions for tem- 
perature variations. In this section we explore the effects ot 
compositional variations parameterized in terms of the vari- 
ation in the Mg # in the upper mantle along the ndge. Our 
motivation for parameterizing compositional variations sim- 
ply in terms of Mg # is that differences in this quantity yield 
significant variations in seismic velocity and density, in 
contrast to most other measures of degree of melt extraction. 

Partial derivatives of density and seismic velocity with 
respect to Mg # are those for olivine, obtained from Akimoto 
[197' ) ] These values were measured on a suite of samples 
ranging from pure forsterite (Mg 2 Si0 4 ) to pure fayalite 
(Fe->Si0 4 ). While these partial derivatives are at standard 
temperature and pressure, it is expected that a change to 
elevated temperature and pressure will have only a second 
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Fig. 14. Same as Figure 13 except for that the viscosity struc- 
ture includes a zone extending from the base of the lid to a depth of 
200 km with a viscosity equal to 0.01 that of the underlying mantle. 
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Fig. 15. Results of combined inversion of geoid, bathymetry, 
and 55-5 travel time residuals for variations in upper mantle 
composition (Mg #). The viscosity structure is taken to consist of a 
40-km-thick high-viscosity lid overlying a constant-viscosity half- 
space. ( a ) Three solutions for along-axis composition variations: 
Dotted line indicates composition perturbations constrained to be 
uniform over 0-150 km depth. Long-dashed line indicates composi- 
tion perturbations constrained to be uniform over 0-300 km depth. 
Short-dashed line indicates composition perturbations constrained 
to be uniform over 0-650 km depth. ( b ) Observed (solid line) and 
predicted along-axis profiles of 55-5 travel time residual, (c) Ob- 
served and predicted along-axis geoid profiles, id) Observed and 
predicted along-axis bathymetry profiles. 
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Fig. 16. Same as Figure 15 except for that the viscosity struc- 
ture includes a zone extending from the base of the lid to a depth of 
200 km with a viscosity equal to 0.0 1 that of the underlying mantle. 


order effect, since temperature and pressure corrections 
work in opposite directions [Jordan, 1979]. Above the soli- 
dus temperature, however, the amount and distribution of 
partial melt, which may depend strongly on composition and 
particularly volatile content, are important. The presence of 
melt is likely to have a larger effect on shear wave velocities 
than on bulk density. Calculations of melt migration, how- 
ever, suggest that once created, melt segregates rapidly by a 
percolation mechanism [e.g., Scott and Stevenson , 1989], so 
that the melt fraction present in the mantle at any given time 


is probably small. Studies of mantle peridotites [Johnson et 
ai , 1990] also support the importance of fractional melting. 

It is straightforward to convert equations (3) and (4) to 
relations between geoid or topography and a compositionally 
induced density perturbation by means of the relation 

A p = -p 0 aA7 (13) 


Compositional anomalies at depth yield a dynamic topogra- 
phy h(k) given by 


AA(*) t P™ H(k, z) AMg(&, z) dz 

P o ~ Pw 5Mg J,^ 

(14) 


where A Mg represents the fractional change in the Mg #. 
Compositional anomalies yield a geoid anomaly 


A N(k) = 


27lT Bp Cz* 
gk 8Mg I 

•* «.mi 


G(k , z) AMg (*, z) dz 


(15) 


For a compositional perturbation at depth the resulting 
two-way travel time perturbation is given by 


A t(k) = 2 


dv s 

a Mg 



AMg ( k , z) dz 


1/2 


(16) 


Using equations (14H16), an inversion scheme similar to 
that used for thermal perturbations is formed. The solution 
vector now has the form 

m = [AMg ■■■, AMg (*„)] r (17) 


The data vector remains the same as in equation (6), while 
the matrix of coefficients, A, changes to reflect the relation 
between the data and mantle composition, rather than tem- 
perature, as outlined in equations (14M16). 

The results of the inversions for compositional variations 
are summarized in Table 2 and in Figures 15 and 16. We are 
unable to match simultaneously both 55-5 travel time 
residuals and geoid and bathymetric anomalies with solely 
mantle compositional variations for either a constant- 
viscosity mantle or one with a low-viscosity zone. This is not 
surprising as the travel times are for the most part positively 
correlated with geoid and bathymetry, but compositional 
variations (at least for the Mg 2 Si 0 4 -Fe 2 Si 0 4 system exam- 
ined here) have an opposite effect on travel time and 
geoid-bathymetry. The xl values for differential travel time 
(Table 2) are extremely large (>4), indicating that the fitting 
function poorly describes the travel time data. 

For the constant-viscosity case, the fit to the geoid is 
excellent, and the fit to bathymetry is slightly better than in 
the inversion for temperature. The fit to 55-5 residuals is so 
poor that the variance reduction is negative for travel time. 
Large compositional changes would be required to affect 
travel times, whereas only small compositional changes are 
needed to produce significant density contrasts to match the 
geoid signal. The total variance reduction and xl values for 
the constant viscosity case do not vary greatly (variance 
reduction from 32 to 33%, xl values from 2.60 to 2.63) for 
compositional changes constrained to be over different depth 
intervals, although the variance reductions and xl values for 
individual data sets (bathymetry, geoid, travel time) vary 


19,998 


Sheehan and Solomon: Joint Inversion of Travel Time, Geoid, and Depth 


CONSTANT VISCOSITY MANTLE 


TEMPERATURE 



LATITUDE. 6 N 



Fig- 17. Results of combined inversion of geoid, bathymetry, and SS-S travel time residuals for both upper mantle 
temperature and composition variations. The viscosity structure is taken to consist of a 40-km-thick high-viscosity lid 
overlying a constant-viscosity half-space, (a) Three solutions for along-axis temperature variations: Dotted line 
indicates composition perturbations constrained to be uniform over 0-150 km depth. Long-dashed line indicates 
composition perturbations constrained to be uniform over 0-300 km depth. Short-dashed line indicates composition 
perturbations constrained to be uniform over 0-650 km depth. ( b ) The corresponding solutions for along-axis 
composition variations, (c) Observed (solid line) and predicted along-axis profiles of SS-S travel time residual, {d) 
Observed and predicted along-axis geoid profiles, (e) Observed and predicted along-axis bathymetry profiles. 


significantly from model to model (see Table 2). The range in 
Mg # is about 1% if the variation is constrained to the depth 
range 0-150 km and only 0.1% for the 0-650 km depth range. 

Figure 16 shows inversion results for the model with a 
low-viscosity zone. A good fit to both geoid and bathymetry 
is found, although the alignment of predicted and observed 
geoid is not as good as in the constant-viscosity case. The fit 
to bathymetry is the best of any models so far. The total 
variance reduction is still low (43—49%) and total xl va l ues 
are high (2.0-2.2), due to the poor fit to travel times (negative 
variance reduction in all cases except the 0-650 km model). 


The range in Mg # is 2.4% if constrained to 0-150 km depth, 
1 .3% over 0-300 km depth, and 0.5% over 0-650 km depth. 

Joint Inversion for Temperature and Composition 

We next explore whether a combination of temperature 
and compositional variations can provide a good match to 
the observed geoid, travel time, and bathymetry. Joint 
inversions provide improved fits to all data at the expense of 
introducing additional free parameters. For these inversions 
the data vector remains the same as in equation (6), the 
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Fig. 18. Same as Figure 17 except for that the viscosity structure includes a zone extending from the base of the lid 
to a depth of 200 km with a viscosity equal to 0.01 that of the underlying mantle. 


solution vector is modified to include both temperature and 
composition, and the matrix of coefficients, A, includes the 
effects of both temperature and composition. The matrix- 
building equations become, for example, for topography, 

AM*) = P ° g -- P”“ H{k, z ) A T{k, z) dz 

Po ~ P” Jzw 

+ — f H(k , z) AMg ( k , z) dz (18) 

P o - Pw dMg J ZmM 

which is simply a combination of equations (3) and (14). The 
new geoid equation comes from a combination of equations 
(4) and (15) and the travel time equation from a combination 
of equations (5) and (16). Cross terms, such as compositional 
changes induced by increases or decreases in temperature, 
are neglected. 


The travel time residuals are perfectly fit in the joint 
inversions for temperature and composition (Table 2). This 
occurs because of the way the model parameters act in a 
similar manner on both geoid and bathymetry, producing a 
singular matrix if only geoid and bathymetry data are in- 
verted for both temperature and composition. Undamped 
least squares always provides perfect solutions when the 
number of equations is equal to the number of unknowns 
unless the matrix to be inverted is singular. If we perform an 
inversion including only travel time and geoid data, we have 
the same number of equations as unknowns, the matrix is 
nonsingular, and we obtain perfect fits to both travel time 
and geoid. Similarly, if we perform an inversion of travel 
time and bathymetry data, we again obtain perfect fits to 
both data sets. If we perform an inversion of geoid and 
bathymetry data, however, we are unable to obtain solutions 
without applying damping. In the joint inversion of travel 
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Fig. 19. Same as Figure 17 but compositional variations are constrained to be from 0-50 km depth only. 


time, geoid, and bathymetry, we have more equations than 
unknowns, and the inversion is overdetermined. However, 
the travel times are perfectly determined in this case because 
of the nonuniqueness inherent with geoid and bathymetry. 
We have performed undamped inversions with various 
weightings on the geoid, bathymetry, and travel time data, 
and in all cases the travel times remain perfectly fit. 

The results for the joint inversion for temperature and 
composition are summarized in Table 2 and Figures 17 and 
18. The geoid residuals are well-modeled in all cases. The 
topography is best fit for the case with a low viscosity zone. 
Resolution of the depth interval of the most important lateral 
variations is rather poor. The topography is fit marginally 
better for the case where temperature and compositional 
anomalies are constrained to be shallower than 300 km. For 
the constant viscosity mantle, the temperature variations 
range from 210 K, if constrained to 0-150 km depth, to 55K 
if over 0-650 km depth; variations in Mg # range from 1.5% 
if over 0-150 km depth to 0.4% for 0-650 km depth. For the 
case with an upper mantle low viscosity zone, the tempera- 
ture variations are similar to those in the constant viscosity 


case, but the variations in Mg # are larger, from over 2% for 
0-150 km depth to nearly 1% for 0-650 km depth. 

We have also performed joint inversion for temperature 
and composition with Mg # variations constrained to be in 
the upper 50 km of the lithosphere so as to mimic composi- 
tional variations due solely to variable melt extraction at the 
ridge. Temperature perturbations were allowed to remain 
within the depth ranges adopted earlier. The results for these 
inversions are summarized in Table 2 and Figures 19 and 20. 
The variance reduction was similar for the constant viscosity 
case and for the model with a low-viscosity zone. In general, 
the geoid is fit very well, the predicted amplitudes are a bit 
low for travel time residuals, and the topography fit is 
slightly out of phase. For the constant viscosity mantle, the 
range in temperature is 210 K over 0-150 km depth and 25 K 
over 0-650 km, while Mg # variations constrained to be 
confined to 0-50 km depth were over 5%. For the case with 
a low-viscosity zone, the temperature variations were not 
dramatically different from those in the constant viscosity 
case, and variations in Mg # were about 4.5%. The inversion 
solution shows high temperatures near 30°N and low tern- 
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Fig. 20. Same as Figure 18 but compositional variations are constrained to be from 0-50 km depth only 


peratures in the region from 50° to 60°N. Iceland also 
appears to be underlain by high-temperature mantle. Going 
from south to north along the ridge, compositional variations 
indicate low Mg # in the vicinity of 20°-30°N , high Mg # in 
the Azores region (40°N), low values near 50°N, and high 
values near 60°N. 

Discussion 

The temperature and compositional variations in Figures 
13-20 are broadly consistent with observed travel time, 
geoid, and bathymetry anomalies in the North Atlantic 
region. Temperature variations alone can account for most 
of the observed anomalies. In contrast, compositional vari- 
ations alone cannot match all anomalies simultaneously. We 
infer that a component of the observed anomalies is due to 
long-wavelength variations in upper mantle temperature. 
Joint inversions for temperature and composition provide 
better fits than single-variable models but at the expense of 
introducing additional free parameters. 

It is difficult to select a “best” model from the suite of 
inversions presented. The variance reductions and xl values 


in Table 2 serve as guides, but independent criteria may 
allow us to reject some of the models, even those with high 
variance reductions and xl values near 1. In particular, 
those models with large temperature variations (well in 
excess of 100 K) can be seriously questioned. Lateral 
temperature variations at upper mantle levels beneath oce- 
anic ridges are thought to be no more than about 300 K 
globally [Klein and Langmuir , 1987], so a variation in 
temperature of 230 K (as in the inversion with a low- 
viscosity zone and a 150-km-thick layer of temperature 
perturbations) solely within a section of the North Atlantic is 
probably unreasonably large. Further, as McKenzie [1984] 
and White and McKenzie [1989] have noted, relatively small 
increases in mantle temperature above values typical for the 
mid-ocean ridge are sufficient to cause large increases in melt 
production. Their models indicate that for fixed bulk com- 
position, an increase of 100 K above normal doubles the 
amount of melt while a 200 K increase can quadruple it. Such 
increased melt production should lead to approximately 
corresponding increases in crustal thickness. Variations in 
oceanic crustal thickness away from fracture zones, how- 
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ever, are generally thought to be small, with thicknesses 
typically 6-7 km and ranging from 4.5 to 8.5 km [Spudich 
and Orcutt , 1980; White , 1984; Purdy and Detrick, 1986]. In 
the joint inversion for temperature and composition, temper- 
ature variations if confined to 150 km depth are excessive 
(over 200 K) and if the variations extend over 0-650 km the 
fit to topography is poor, especially for the constant- 
viscosity mantle. On the basis of these results we prefer the 
models with temperature variations occurring over 0-300 km 
depth. For the constant viscosity mantle, the temperature 
variation is 110 K, and the variation in Mg # is 0.75%. For 
the case with an upper mantle low-viscosity zone, the 
predicted temperature variation is 125 K, and the variation in 
Mg # is 1.1%. The total variance reduction is greater in the 
model with a low-viscosity zone. 

Even a temperature variation of about 100 K is high for a 
mantle of constant composition, since we do not observe 
increased crustal thickness in regions that our models indi- 
cate have high temperatures. The assumption of approxi- 
mately constant upper mantle composition warrants discus- 
sion. In particular, lateral variation in trace amounts of 
mantle volatiles may have a large effect on seismic velocity 
at a given temperature. The presence of even a slight amount 
of water, for instance, is sufficient to cause a significant 
decrease in the initial melting temperature of peridotite 
[Wyllie, 1971]. Estimates of volatile contents and their lateral 
variations in the North Atlantic region have been made from 
measurements of abundances of halogens, Si0 2 , K : 0, and 
H 2 0 in basalts and from the volumes of vesicles in basalts 
[Schilling et al . , 1980, 1983; Schilling , 1986; Michael , 1988]. 
These studies indicate that Cl, Br, F, and H : 0 contents 
increase toward the Azores and Iceland and that H : 0 is 2-3 
times more abundant in Mid-Atlantic Ridge basalts erupted 
over the Azores platform than at adjacent normal ridge 
segments. The effect of volatiles on density and shear wave 
velocity will be slight at subsolidus temperatures but can be 
major over the melting interval [ Goetze , 1977]. The presence 
of melt will act to decrease significantly the seismic velocity 
[Duschenes and Solomon, 1977] and, to a lesser extent, 
lower the density of the mantle. To the extent that seismic 
velocity depends on proximity of the temperature to the 
solidus temperature [Sato et al., 1988, 1989], volatile content 
can trade off with temperature in its effect on velocity at 
subsolidus conditions. Thus, variation in volatile content 
could lessen the variations in melt production implied by the 
inversion solutions. 

Even without significant variations in volatile content, it is 
clearly an oversimplification to parameterize mantle compo- 
sition in terms of only a single quantity. Further, we have 
assumed that the partial derivatives of bulk density and 
seismic velocity with respect to Mg # that are those for 
olivine [Akimoto, 1972]. The work of Jordan [1979] indicates 
that these derivatives remain nearly constant for many 
different mantle compositions (i.e., pyrolite-type composi- 
tions with various amounts of olivine, orthopyroxene, cli- 
nopyroxene, spinel, and garnet), so the latter assumption is 
sound. However, at any given Mg #, orthopyroxene and 
clinopyroxene have lower velocities and are less dense than 
olivine, while garnet and spinel are seismically faster and 
denser than olivine [ Jordan , 1979], so an increase in the 
weight percent of orthopyroxene and clinopyroxene or a 
decrease in the weight percent of garnet and spinel with 
respect to olivine in the mantle could counteract some of the 


temperature variations obtained under the assumption of 
effectively uniform, mineralogy. Several studies [Wood, 
1979; Jaques and Green , 1980; Dick et al., 1984] have 
suggested that compositional variations in the mantle are 
plausible. Indeed a number of workers [e.g., Davies , 1984; 
Allegre et al., 1984] favor dynamic models for the mantle in 
which dispersed heterogeneities of various sizes and shapes 
are passively embedded in a continually mixed, convecting 
mantle. Variations in modal fractions of olivine, orthopyrox- 
ene, and clinopyroxene in peridotites recovered along the 
Mid-Atlantic Ridge have been reported in several studies 
[Dick et al., 1984; Michael and Bonatti , 1985]. These varia- 
tions are typically attributed to different degrees of melt 
extraction but could also be partially due to intrinsic upper 
mantle heterogeneity. For example, the relative fractions of 
olivine, clinopyroxene, and orthopyroxene indicated by 
Michael and Bonatti [1985] at 26°N and 30°N, if extended to 
depth, could counteract a portion of the temperature differ- 
ences indicated by the inversion solutions for these regions. 

Chemical analysis of dredged peridotites in the north 
Atlantic indicate a range of about 2.5% variation in Mg # 
[Michael and Bonatti , 1985]. This value is intermediate 
between what we find for models with compositional varia- 
tions constrained to be shallow (4.5-6% variation) and those 
models with compositional variations in the same depth 
ranges as the thermal variations (1-2%). This suggests tha? 
compositional variations may be concentrated slightly shal- 
lower than the temperature variations. Michael and Bonatti 
[1985] present an along-axis profile of Mg # variations from 
dredged peridotites which can be compared with our calcu- 
lated profile. The main feature in their profile is a zone of 
high values of Mg # in the Azores region, from 34° to 45°N, 
relative to the rest of the ridge, consistent with our modelling 
results. Their data sampling is too sparse to delineate other 
long-wavelength features. Their average value for 26 D N also 
has a high Mg # relative to adjacent data. This is consiste r 
with our observation of early SS-S travel times and Io\> 
geoid in this region. This anomaly is of too short a wave- 
length (<1000 km), however, to resolve in our inversions. 
We should note that such comparisons merit caution, as 
small-scale features, such as those due to ridge segmenta- 
tion, can produce large differences in composition between 
peridotites over scales of tens of kilometers. In addition, 
dredged peridotites are mostly from fracture zone environ- 
ments, which may not be representative of typical ridge 
mantle [Dick, 1989]. 

On the SS-S residual profile the Iceland region appear^ 
a local maximum (late SS) but the Azores hotspot does not 
show a distinct seismic signal. The inversion results for these 
two regions are also markedly different. The results of the 
joint inversion for temperature and composition predict a 
high Mg # in the Azores region while indicated temperatures 
are not anomalously high. At Iceland, in contrast, high 
temperatures dominate. Work by Schilling [1986] and Bon - 
aiti [1990] outlines the differences in geochemical signatures 
between the Azores and Iceland hotspots. These wori; ar s 
suggest that Iceland is a ‘‘traditional”' plume hotspot, v. ... a 
predominantly thermal origin, but that the Azores might be 
more aptly named a “wet spot" because of the presence of 
excess hydrous phases and the lack of a thermal anomaly- 
Langmuir and Hanson [1980] noted low Fe contents for 
basalts from the Azores region and suggested that, in con- 
trast to Iceland, the Azores overlie a mantle that is depleted 
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in basaltic constituents and has been modified by a volatile- 
rich component. Bonatti [1990] suggests that because the 
Azores hotspot is rich in volatiles, enhanced melt production 
could occur with little or no increase in temperature. The 
high Mg # indicated in our inversions allows the region to be 
seismically fast (as we observe) but of low density (as geoid 
and bathymetry require). The results are consistent with the 
hypothesis that the Azores hotspot is not associated with a 
significant plumelike thermal anomaly. Inversion of surface 
wave dispersion data can potentially provide further tests of 
these ideas, but studies to date have yielded apparently 
conflicting results. Results of several such investigations 
[Nakanishi and Anderson, 1984; Tanimoto, 1990; Zhang and 
Tanimoto , 1990, 1991] suggest that the Azores region is 
seismically slow at depths less than 300 km, but a study 
utilizing 50- to 200-s-period Rayleigh waves by Mocquet et 
al. [1989] does not. These differences may be partially 
attributable to the differences in wave periods employed and 
mode of analysis from study to study. It may be possible that 
what appear to be low velocities at the Azores are a result of 
horizontal smoothing of the low velocities along the ridge 
and have little to do with the actual structure in the Azores 
region. Few of these long-period surface wave studies re- 
solve a distinctive anomaly at Iceland. Clearly, more work is 
needed to resolve the upper mantle velocity structure of 
hotspot regions. 

Zhang and Tanimoto [1991] have constructed a profile of 
Love wave phase velocity variation along the Mid-Atlantic 
Ridge at a period of 100 s. Their profile has many features in 
common with our along-axis profile of SS-S differential 
travel time, but we also note several differences. In both 
profiles the region from 10° to 40°N is “slow” relative to the 
region 40°-70°N. The ridge section from 45° to 55 3 N is 
relatively “fast” in both profiles, with velocities becoming 
gradually lower northward to the Iceland platform and 
gradually higher again north of Iceland (65°-75 a N). The most 
notable difference is at 40°N, in the vicinity near the Azores 
hotspot. The model of Zhang and Tanimoto [1991] has low 
velocities (-0.8% less than the average Love wave phase 
velocity), while our SS-S travel time residuals indicate no 
distinct anomaly in this region. This discrepancy could be 
caused by a number of different factors. The surface waves 
and body waves both average over a portion of the upper 
mantle, but each does so in a different manner. The 100-s 
surface waves have kernels which peak at shallow depth 
(0-100 km) and have decreased sensitivity to depths as great 
as about 400 km. The SS waves integrate anomalies through 
the entire upper mantle, with equal resolving power at all 
depths. The lateral averaging of surface wave and body wave 
data are also quite different. If a mantle plume underlying a 
hot spot has a shallow, extended “head” of anomalously hot 
mantle material [e.g., White and McKenzie , 1989], then a 
surface wave passing through the plume center may display 
a greater velocity anomaly than a body wave sampling the 
plume at nearly vertical incidence. By constructing a profile 
from data with bounce points distributed in age from 0 to 100 
Ma, we perform lateral averaging which may further smooth 
out any hotspot signal. However, a map view of the individ- 
ual travel time residuals in the immediate vicinity of the 
Azores (Figure 6) does not show any evidence for a signifi- 
cant anomaly. The Azores region is not anomalously slow in 
the Love wave phase velocity data of Zhang and Tanimoto 
[1991] at periods of 150 s and greater, indicating that the 


source of the anomaly in the 100-s phase velocity most likely 
arises from the shallow mantle. 

Bonatti [1990] has constructed profiles of the equilibrium 
temperature of dredged peridotites along the Mid-Atlantic 
Ridge axis from 0 to 60°N by means of two different 
geothermometers [Wells, 1977; Lindsley, 1983]. Comparison 
of these profiles with the along-axis temperature variations 
obtained from our inversions reveals a number of qualitative 
correlations as well as a few discrepancies. The range of 
temperature variations in the profile based on the Lindsley 
[1983] geothermometer is about 150 K, neglecting high 
values termed “anomalous.” When the high values are 
included the range increases to 400 K. The profile utilizing 
the Wells [1977] geothermometer has a range of 100 K 
neglecting the anomalous values and 350 K including them. 
The highest temperatures in our inversions are near 30°N 
(Figures 17 and 18), a region showing a slight peak in 
Bonatti’s temperature profile estimated according to Linds- 
ley [1983] and a very weak rise in the profile utilizing the 
Wells [1977] geothermometer. There is a small dip in tem- 
perature at 26°N (a region which we find to be seismically 
fast) in the Lindsley [1983] and Wells [1977] profiles, but the 
difference may not be significant considering the errors. 
Bergman and Solomon [1989] also found the upper mantle 
near 26°N to be seismically fast from an analysis of teleseis- 
mic P wave travel time residuals from earthquakes in this 
region recorded by a local ocean bottom seismic network. 
The lowest temperatures on the profiles of Bonatti [1990] are 
at 43°N. Temperatures from our inversion solutions are also 
low in this region, although the Bonatti [1990] profiles 
indicate an increase in temperature proceeding north from 
43°N to 53°N, whereas our results favor continued low 
temperatures. Part of the difference between our results and 
the geochemical studies may be attributed to the fact that the 
depth sampled by basalts and peridotites is likely to be 
shallower than the layer thicknesses of most of our models. 
Assuming that the 6-km-thick oceanic crust was formed by 
9-22% partial melting of the mantle [Klein and Langmuir , 
1987], then the volume of residual peridotite will extend from 
the base of the crust to somewhere between 30 and 70 km 
depth. The amount of depletion will vary with depth if we 
assume a fractional melting model. Our models with compo- 
sitional variations confined to depths less than 50 km are 
most representative of shallow fractionation and differenti- 
ation. 

Alternate petrological measures of upper mantle temper- 
ature may be obtained from the major element chemistry of 
mantle-derived basalts. Klein and Langmuir [1987, 1989] 
suggest that the quantity Fe 8 0 (basalt FeO content corrected 
to 8.0 wt % MgO) is related to mantle temperature in the 
source region, with low Fe 8 0 indicating a relatively cool 
mantle and a narrow depth interval of melt generation and 
high Fe g>0 indicating a hotter mantle and a greater depth 
interval of melt generation. A comparison of the temperature 
variation obtained in our inversions with an along-axis 
profile of Fe 8 o [Klein and Langmuir , 1989] shows several 
striking correlations, particularly peaks in both quantities 
near 30°N and near Iceland. 

Several improvements in future studies of the type pre- 
sented here may be envisioned. Our models thus far have 
been limited to simply parameterized one-dimensional vari- 
ations in temperature and composition within a single layer. 
It is likely that these lateral variations are not constant 
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within a given layer and that there are two-dimensional 
lateral variations independent of lithospheric aging. The 
techniques outlined in this paper can be generalized to a 
multilayer system and to two-dimensional wavenumber (see 
equation (8)), but we do not feel that the resolution of our 
data can justify more complicated models at this time. 
Kernels for seismic surface waves are strongly peaked in the 
upper mantle, and such data would provide a useful con- 
straint in future models. The inclusion of surface wave data 
would help to distinguish between lithospheric and astheno- 
spheric effects and may allow for two or more independently 
resolved layers. Extension of the modelling to three dimen- 
sions would permit an assessment of the degree to which 
mantle anomalies beneath the ridge extend off axis. Implicit 
in our age correction is the assumption that the anomalous 
properties of the ridge mantle are steady state on a time scale 
of 100 Ma. Recent seafloor surveys and theoretical studies 
[e.g., Pockalny et al. , 1988; Scot! and Stevenson, 1989] bring 
this assumption into question and suggest that at least on 
short time scales (<1 m.y.) and at slow spreading rates (as in 
the Atlantic) intermittent periods of melting and crustal 
formation may be separated by periods with little or no melt 
production. These temporal variations are likely to be aver- 
aged out, however, over the typical horizontal wavelength 
(> 1 00 km) of a long-period SS wave. 

Another limitation of our models is that they depend on 
the assumed values of several physical constants. It is 
straightforward, however, to estimate the effect of choosing 
different values. The viscosity structures that we employ are 
also quite simple but have been chosen to represent two 
models widely invoked in other studies, a constant or nearly 
constant viscosity mantle [e.g., Peltier , 1989] and a mantle 
with a thin low-viscosity layer [e.g., Craig and McKenzie , 
1986; Robinson et al ., 1987]. The viscosity structure of the 
Earth may be temperature and pressure dependent or vary 
laterally, but we have not considered viscosity structures of 
this type [e.g., Revenaugh and Parsons, 1987]. Much work 
remains to be done to determine ways to incorporate lateral 
variations in viscosity into quantitative treatments of these 
problems. We have not modelled the effects of partial 
melting which could accompany the temperature variations 
that we predict. The effect of retained melt on the physical 
properties of the mantle depends critically on the melt 
fraction and geometry, characteristics presently poorly 
known. Sato et al. [1988, 1989] downplay the importance of 
partial melt and suggest that most mantle seismic velocity 
anomalies can be explained by temperature variations at 
subsolidus conditions. The combined analysis of both shear 
and compressional differential travel times also suggest that 
significant partial melting is not required to explain the travel 
time residuals in the North Atlantic region [Woodward and 
Masters , 1991a]. 

Conclusions 

We have measured nearly 500 55-5 differential travel 
times for paths in the North Atlantic region. The 55-5 travel 
time residual decreases linearly with square root of age, in 
general agreement with the plate cooling model, to an age of 
80-100 Ma [Parsons and Sclater , 1977]. Azimuthal anisot- 
ropy is not clearly resolved; in particular, the azimuthal 
patterns of our data are not consistent with the preferred 
upper mantle anisotropy model of Kuo et al. [1987] for the 


North Atlantic. If azimuthal anisotropy is present in the 
North Atlantic, it is not uniform over the entire region. An 
along-axis profile of age-corrected travel time residuals 
displays significant variations at wavelengths of 1000-7000 
km. The largest of these variations are robust with respect to 
selective removal of portions of the data. 

We have formulated a joint inversion of travel time 
residuals and geoid and bathymetric anomalies for lateral 
variation in upper mantle temperature and composition. On 
the basis of variance reduction, inversion for temperature 
favors the presence of an upper mantle low-viscosity zone 
and temperature anomalies concentrated at depths less than 
300 km. We are unable to match travel time residuals 
simultaneously with geoid and bathymetry solely with lateral 
variations in bulk composition (Mg #). Joint inversions for 
temperature and composition provide good fits to both travel 
time and geoid regardless of viscosity structure or layer 
depth and thickness, but the best fits to bathymetry come 
from models with a low-viscosity zone and thermal or 
compositional variations confined to shallow depth. The 
Mg # variations predicted in the joint inversion for temper- 
ature and composition are comparable to those found by 
Michael and Bonarti [1985] in a study of dredged peridotites 
along the Mid-Atlantic Ridge and may be related to varia- 
tions in melt production along the ridge. 

The preferred inversion solutions have variations in upper 
mantle temperature along the Mid-Atlantic Ridge of about 
100 K. For a constant bulk composition, such a temperature 
variation would produce about a 7-km variation in crustal 
thickness [White and McKenzie , 1989], larger than is gener- 
ally observed [Spudich and Orcutt , 1980; White , 1984; Purdy 
and Detrick , 1986]. Introducing compositional variations as 
well as temperature variations in the inversions does not 
change the range of temperature appreciably. The presence 
of volatiles in the mantle can have a strong effect on 
temperatures required for melting, and variations in volatile 
content along the ridge may reduce the large variation in 
melt production implied by the lateral temperature variations 
indicated in our models. 

Appendix A: Estimation of Errors for 55-5 
Differential Travel Times 

It is important to quantify the uncertainties in the differ- 
ential travel time measurements. After cross-correlation, the 
“quality” of each individual 55-5 measurement is rated and 
a grade is assigned. The cross-correlation coefficient, which 
describes the degree of fit between the synthetic and real 55 
phases, is used as an objective aid in the assignment of 
quality. However, our final assignment of quality is largely 
subjective and based upon visual inspection of the “synthet- 
ic” 55, real 55, and cross correlogram, taking into account 
the sharpness of the arrivals and their alignment, the clarity 
of the seismogram, and the appearance of a single clear peak 
in the cross-correlation function. An “A” quality grade 
indicates an excellent fit, “B” quality indicates good phase 
alignment but only a fair fit, and a “C" quality grade 
indicates a poor fit or some ambiguity as to phase alignment. 
In addition to A, B, and C grades, there were data that were 
rejected due to poor signal to noise ratio for either the 5 or 
55 phases. 

Assuming that the uncertainty in an individual measure- 
ment comes from a combination of measurement error, 



Sheehan and Solomon: Joint Inversion of 


Travel Time. Geoid. and Depth 


:o,oo5 


unmodeled lower mantle structure, and epicentral error, we 
write, for example, for the measurement variance of an “A” 
quality datum: 

°A = *Am + O-lm + ^Pi (A1) 

vhere cr A is the total uncertainty, <r Am is the measurement 
error, cr\ m is the uncertainty due to unmodeled lower mantle 
structure"! and c r epi is the epicentral error. We assume that 
CT| m and o- epi are the same for A, B, and C quality measure- 
ments, but the measurement error is obviously a strong 
function of data quality. 


Effect of Epicentral Error 

In general, epicentral errors affect the travel times only 
slightly. The events used in this study were well recorded by 
a large number of stations over a wide range of azimuths, 
and typical epicentral mislocations are probably less than 10 
km (which would yield a differential travel time error of 0.35 
s at 75° distance). The travel times are even less sensitive to 
errors in focal depth; an error in depth of 25 km contributes 
only about 0.3 s to the 55-5 residual. Using the rule of 
thumb that one standard deviation is about one half of the 
estimated extremes, we adopt cr epi = 0.75 s as a conservative 
estimate of epicentral error. 


Effect of Unmodelled Lower Mantle Heterogeneity 

We estimate the likely magnitude of lateral variations in 
the shear wave velocity of the lower mantle from models of 
lower mantle heterogeneity in P wave velocity (such as 
model L02.56 of Dziewonski [1984]). The average variation 
in travel times of direct P waves bottoming in the lower 
mantle is in the range ±0.5 s. Global tomographic studies by 
Dziewonski and Woodhouse [1987] indicate that the scaling 
ratio ( 5 V s IV s )l( 8 V p !V p ) ~ 2 in the lower mantle. Such a 
scaling is also suggested by comparison of lower mantle P 
wave models with the recent lower mantle 5 model of 
Tanimoto [1990]. Assuming such an 5 to P velocity anomaly 
scaling, the resulting variation in 5 wave arrival time con- 
tributed by the lower mantle would likely be about ± 1.5 s, a 
fraction of the observed range in 55-5 residual. While the 
major features of lower mantle model L02.56 [Dziewonski, 
1984] and the lower mantle portions of Tanimoto s [1990] 
model are for the most part similar, enough differences exist 
that the application of a lower mantle “correction” to our 
data might add more uncertainty than it removes. Further, 
absolute S wave travel times do not show enough variance 
for us to suspect large lower mantle effects [e.g., Randall, 
1971; Girardin and Poupinet, 1974; Hart and Butler, 1978; 
Uhrhammer, 1978, 1979], and the work of Gudmundson et 
al. [1990] indicates that most of the variance from the ISC 
tables is attributable to the shallow mantle, i.e., most of the 
Earth’s heterogeneity is in the upper mantle and the lower 
mantle is fairly homogeneous. On the basis of the above 
information, we set <r [m = 0.5 s for our study. We note, 
however, that knowledge of global-scale structure in the 
lower mantle is still incomplete. Some recent studies [Wood- 
ward and Masters, 1991a, b] suggest that a larger value ot 
cr, is warranted. However, the fact that our study involves 
only a limited geographic area may confine the lower mantle 
variation to only a fraction of the full global range. 


Measurement Error 

As an objective means to obtain error estimates, we 
examine the scatter in A, B, and C quality picks in a small 
region. We measured the root mean square (rms) difference 
between travel time residuals of the same grade (A, B, or C) 
with bounce points separated by less than 80 km and with 
differences in path azimuth at the bounce point of less than 
10° An 80-km distance is less than the horizontal wavelength 
of 55 (which is about 180 km at 25-s period) so we do not 
expect much contribution to the rms difference from actual 
lateral variations in structure. The rms difference for the 16 
\ quality residual pairs which were within 80 km of each 
other was 1.15 s. For B quality picks, an rms difference of 
2 08 s was measured using 20 residual pairs, and for C 
quality picks, 44 residual pairs yielded an rms difference of 
2.96 s. 

We interpret these estimates of the rms differences as 
representing the average overall errors in the A, B, and C 
grade measurements. (Unmodelled lower mantle structure 
should be nearly identical for data with bounce points within 
80 km and at similar azimuths.) Under this interpretation we 
can write, for A quality residuals 

O-irms = O'Am ^ (A2) 

Substituting values of o- Arms and cr ep , into (A2) yields cta* - 
0.87 s. Similarly, for B and C quality measurements, we hnd 
u = 1.94 s and cx Cm = 2.86 s. From (Al), the total 
uncertainty for A, B, and C quality is, respectively, <r A = 
1.25 s, o- B = 2.14 s, and <r c = 3.00 s. 

In the weighted regression experiments the A, B, and C 
quality measurements are weighted inversely by their mea- 
surement variance. 


Appendix B: Errors in the Along-Axis Profiles 

and Construction of Data Covariance Matrix 

Errors in Bathymetry, Geoid, and Travel Time Profiles 

Uncertainties in the along-axis profiles of geoid, bathym- 
etry, and travel times are important information in the 
inversion. The gridded bathymetric data [U.5. Naval Ocean- 
ographic Office, 1985] include corrections for the deviation 
of water column acoustic velocity from the assumed value of 
1500 m s' 1 . The geoid data, provided in the form of a0.25 
x 0 25° grid [Marsh et al., 1986], include corrections for 
orbit errors, instrument and atmospheric propagation ef- 
fects, and solid Earth and ocean tides. 

We have averaged the bathymetry and geoid height values 
within a 1° x 1° box centered at each SS bounce point. The 
averaging yields a representative value for a region over 
approximately one horizontal wavelength of the SS wave 
and acts to smooth out short-wavelength variations. Both 
bathymetry and geoid are corrected for subsidence with 
seafloor age, using the plate cooling model [Parsons and 
Sclater , 1977; Parsons and Richter , 1980]. Error introduced 
into depth and geoid anomalies by isochron mislocation is 
difficult to estimate precisely, but for an error in age of 2 
m.y., depth and geoid errors at 80 Ma would be about 30 m 
and 0.2 m, respectively, while at 2 Ma, an error in age of 2 
m.y. would have a much larger effect, giving depth and geoid 
errors of 350 and 0.3 m. The magnitude of this error 
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highlights the importance of accurate age determination, 
especially at young ages. 

The presence of oceanic sediments is another source of 
error. In the Atlantic Ocean, the sediment thickness in- 
creases regularly from less than 100 m along the Mid- 
Atlantic Ridge toward continental margins where it can 
exceed 1 km [Ewing et al. , 1973; Tucholke , 1986]. A l-km 
sediment thickness leads to corrections to residual depth and 
geoid of about 500 and 0.3 m, respectively [Cazenave et al ., 
1988; Sheehan and McNutt. 1989]. On Atlantic lithosphere 
of 100 Ma age or less the sediment thickness is less than 500 
m in most areas. Hence, neglecting the sediment loading 
correction should not be crucial in this region. 

The along-axis profile of SS-S residual (Figure 10) is a 
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8-14 km thick [ Bjornsson , 1983]. By simple isostatic mass 
balance, the depth anomaly due to excess crustal thickness 
in the Azores region would be about 400 m, and at Iceland, 
200 m to 1.6 km. In general simple variations in crustal 
thickness are insufficient to produce a significant SS-S 
residual. For crustal and mantle S wave velocities of 3.5 and 
4.4 km s'\ a 2-km variation in crustal thickness would 
contribute less than 0.2 s to an SS-S differential travel time 
corrected for differences in bathymetry. However, at Ice- 
land, where the crust is estimated to be as much as 14 km 
thick, the additional SS-S travel time could be up to 0.8 s. 

Data Covariance Matrix 
The data covariance matrix R dd is of the form 


0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 0 

0 0 0 0 

0 cr N (k„) 2 0 0 0 

0 0 < r ,( fc ,) 2 0 0 

00 0 ••• 0 

0 0 0 0 cr,(k„) 2 



weighted moving average of 10 adjacent data points grouped 
by latitude, using the weights discussed in Appendix A. The 
same weights and moving average are applied to geoid and 
bathymetry values at a given SS bounce point, even though 
bathymetry and geoid data sets are each presumed to be of 
uniform quality, in order that these profiles will be consistent 
with the SS-S residuals. The standard error of the mean 
values for SS-S residual range from 0.2 to 1.6 s. For 
bathymetry the range of standard deviations from the mean 
value is from 24 to 370 m and for geoid is 0.08 to 1 .0 m. The 
largest variances in the bathymetry and geoid data come 
from the Iceland region (north of 60°N), and may be due to 
the more complicated tectonics of this region [White, 1988]. 

Before Fourier transforming, the along-axis profiles must 
be interpolated to a constant spacing. We use a simple linear 
interpolation scheme to estimate values at 0.5° spacing. We 
estimate that the typical error in the interpolated data is 
comparable to that in the along-axis moving averages, which 
for bathymetry is of the order of 125 m, for geoid 0.4 m, and 
for travel time 1 s. 


Effect of Crustal Thickness Variations 

Our poor fit to topography in the inversion experiments 
can be at least partially attributed to unmodelled effects such 
as crustal thickness differences. Variations in oceanic crustal 
thickness about the typical value of 6-7 km [Spudich and 
Orcutt , 1980; White, 1984; Purdy and Detrick , 1986] are 
generally thought to be small at horizontal scales of 100 km 
and greater. However, the crust beneath the Azores plateau 
is estimated to be between 8 and 9 km thick [ Searie , 1976; 
Whitmarsh et al. , 1982] and that beneath Iceland is at least 


where erg, erh, and cr, 2 are the nominal variances of the 
bathymetry, geoid, and travel time data, respectively. We 
may choose to construct the data covariance matrix not to 
reflect the true variance of the data but rather to allow 
weighting between the different data sets. In this way, the 
data covariance matrix can be altered to test the relative 
contributions of different data sets to the inversion results. 

In all of our inversions, the covariance matrix is con- 
structed to weight the three data sets approximately equally. 
For example, examination of Figure 10 indicates that at 
3000-km wavelength the amplitude of the geoid signal is 
approximately 4 m, bathymetry 1 km, and travel time 2 s. 
Thus if a value of 1 m is chosen for then a value of 0.5 
s for a t and 0.25 km for o h should yield approximately equal 
weighting of data sets. The corresponding 1/cr 2 values are 
then 1 for geoid, 4 for travel time, and 16 for bathymetry. 
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